{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths for loading/storing results.\n",
    "name = 'results/MixedGaussian_MINEE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the mixed gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.mix_gaussian import MixedGaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 400   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "\n",
    "rep = 1             # number of repeated runs\n",
    "d = 1               # number of dimensions for X (and Y)\n",
    "\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "mg = MixedGaussian(sample_size=sample_size,rho1=rho,rho2=-rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = mg.data\n",
    "        X[i,:,j] = data[:,0]\n",
    "        Y[i,:,j] = data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X/wZXV93/HnG5AN8iNKoCEoXzBoTCzDqN36Y2qijSRCJElDY2riiK0Td4hJJ2bsahXiT8wodBgz0TpuK9E1/mzF6AgGsG1CbRbr4hBcRFvWuAvCgrqSBaTg1333j3s+y/mePb9/3M85574eMwx7v/fccz7n3Hs+7/P5be6OiIjIEbETICIi46CAICIigAKCiIgkFBBERARQQBARkYQCgoiIAAoIEoGZ/bWZ/e6SjvV7ZnaPmT1gZj9RY/tvmdk5y0jbGJjZB83s0tjpkHFQQJBBJBnrQ0lGfI+Z/bmZHddwH2eYmZvZUS3T8BjgCuCX3f04d/9em/2U7N/N7Ml97lMkJgUEGdKvuvtxwDOBfwpcsuTj/yTwY8CtSz6uyCQpIMjg3P3bwOeBs7LvmdkRZnaJme0xs3vNbLuZ/Xjy9g3J/+9LShrPzfn8JjN7t5ndlfz37uRvPwN8I/X5/56XNjN7eXLs75nZxZn3nmVmO8zsPjO728zeY2ZHJ++FtP1dkrZ/ZWaPN7PPmdl3zOz7yb+fWHRdzOz1ZvZtM7vfzL5hZi+sOm7yvpvZq83s/yaffbuZnZl85oCZfTKVzheY2Z1m9kYz+25ScntZSZrON7Obk2P/rZmdXZVemRF313/6r/f/gG8B5yT/Po3FU/rbk9d/Dfxu8u9XArcDPw0cB1wFfDh57wzAgaNKjvM24EbgHwEnA3+bOk7p54GnAQ8AvwBsYlG9tJ5K9z8BngMclezrNuA1qc878OTU658A/iXwWOB44L8Af1lw7KcCdwCnptJ6ZoPjfhY4AfjHwMPAf0uu4Y8DXwNekWz7guScrkjO8fnAg8BTk/c/CFya/PuZwL3As4EjgVck3+OmsvTqv/n8pxKCDOkvzew+4IvA3wB/krPNy4Ar3P2b7v4A8AbgpQ3aDV4GvM3d73X37wBvBV5e87O/CXzO3W9w94eBPwYOhjfd/SZ3v9Hd1939W8D7WWSoudz9e+7+KXf/gbvfD7yjZPsfschon2Zmj3H3b7n77gbHfZe7H3D3W4FdwHXJNfwHFqWxZ2S2/2N3f9jd/wa4GvitnDS9Cni/u3/J3X/k7h9iEWyeU5ZemQ8FBBnSv3D3x7n76e7+and/KGebU4E9qdd7WDwZ/2TNY+R9/tQGn70jvHD3B4FDDc9m9jNJtc8+MzvAIqCdVLQzM3usmb0/qYI6wKLK63FmdmR2W3e/HXgN8BbgXjP7uJmd2uC496T+/VDO63QD/veTcwuKrtHpwGuT6qL7kmB+GotSQWF6ZT4UECS2u1hkRMEaiyqOe1hUjbT5/F01j303iwwPWGToLKp9gvcBXwee4u4nAG8ErGR/r2VRtfLsZPtfCLvO29jdP+ruz0vS78C7Wh63yuPN7NjU66JrdAfwjiSIh/8e6+4fq0ivzIQCgsT2MeCPzOxJSbfUPwE+4e7rwHdYVOH8dMXnLzGzk83sJOBNwF/UPPZ/Bc43s+cljbBvY+M9cTxwAHjAzH4W+L3M5+/JpO14Fk/n95nZicCbiw5sZk81s180s03A/0s+96Oax23jrWZ2tJn9PHA+i/aNrP8EXGRmz7aFY83sxWZ2fEV6ZSYUECS2K4EPs6he+XsWmc2/BXD3H7Coh/9fSRXGc3I+fymwE7gF+CrwleRvlZL6998HPsqitPB94M7UJv8O+B3gfhaZ5Scyu3gL8KEkbb8FvBs4Bvgui4buvyo5/Cbgncm2+1g0ir+x5nGb2sfi3O4CPgJc5O5fz27k7jtZtCO8J9n+duBf10ivzIS5a4EckbkysxcAf+Huhd1fRQKVEEREBIgYEMzsx8zsf5vZ35nZrWb21lhpERGRiFVGZmbAse7+gC3mnPki8IfufmOUBImIrLhWk4b1wReR6IHk5WOS/9SgISISSbSAAJAM2LkJeDLwXnf/Utn2J510kp9xxhnLSJqIyGzcdNNN33X3k6u2ixoQ3P1HwNPN7HHAp83sLHffld7GzLYAWwDW1tbYuXNnhJSKiEyXme2p3mokvYzc/T4WE56dm/PeNnff7O6bTz65MsCJiEhLMXsZnZyUDDCzY4BzWAzXFxGRCGJWGf0Ui1GeR7IITJ90989FTI+IyEqL2cvoFg6foldERCIZRRuCiIjEp4AgIiKAAoKIyOhs3b6Drdt3LP24CggiIgJEHpgmIiKPCqWCW/bs3/D68gufu5Tjq4QgIhJJrKqhIioh1LTsSC0iqyfkL7HyGwUEEZEli101VEQBocJYvzgRma9Y+YsCgojIksWuGiqigFBhrF+ciEjfFBBERCIZ2wOmAkJNY/viRORRKsH3Q+MQREQEUAlBRCZMvQD7pRKCiIgAKiGIyISpF2C/VEIQERFAJQQRmQGVDPqhEoKIiAAKCCIiklBAEBERQAFBREQSCggiIgIoIIiISEIBQUSiG9vawqtKAUFERAANTBORiDQ53biohCAiIkDEEoKZnQZsB04BDgLb3P1PY6VHRJZPk9ONS8wqo3Xgte7+FTM7HrjJzK53969FTJOISFQxg2O0gODudwN3J/++38xuA54ADBIQ5voEMtfzktWi3+84jKJR2czOAJ4BfCnnvS3AFoC1tbWlpks2UvCRuYv5Gx9DA3v0gGBmxwGfAl7j7gey77v7NmAbwObNm73p/sdwkYcw1/MSkXiiBgQzewyLYPARd78qZlqk2JDBR4FMxmAMD1hjaGCP2cvIgA8At7n7FUMdZwwXeQjLOq+t23ewe98BzjzlhEH2LyLjEbOE8M+AlwNfNbObk7+90d2viZgmyXHmKSdw+YXPHaRkoCovGUqT39SYHhxjHjtmL6MvAras4801oxmyZAAbM2yVFETmLXqjsgxjiCedUFLow5ieyGReupQ+V/13qIAguZRhy1Tt3ndYZ0WpSQFhILEy0qnVzY81XTINeb/vbLWmfmP1KSCsmKYBQjeTjEnZ7zf7MHTspu7Z29gfqPqmgNCz2E/oquqRZRpTSTjb6UEdIJpTQFgRsQOVSBd1nv777B69qveLAkLPxvKEHuuJbe43jCzEzjDVTjAMBYQVsYxAtXvfAbZu36Gbc+LGGNyLfr956zD3ke6xPNgtmwLCQFblB5R9UlRQWA3LzDDLjqHfWb8UEFbMUCWD4MGH1xUUJip2NVAd2bRUpbHrOYzp3JdBAUE6CY14u/cd4MGH1wH17lglyygZjDlAzY0CgnSWDgoKBtPVZzXQMmbhheJgoWDSjgKC9CLd3U+kD6vasBuTAoL0Tk9l0zaF/vtVwULBpB0FBBEZTB8TzSkzXx5zb7xMcTSbN2/2nTt3xk6GVIjxVKYnwXG54LJrgUc7GOh7icvMbnL3zVXbqYQgIp1l249Cj7OqqiMF8nFRQJDexSgZqN2if02vZageUk+z6VJAkF7VmZ5YmfVwhrrGVU/42Unnzj79xA3bVX2uS7r1u+qPAoJMmnqT9K/petp5I9VBJYUpUqOy9CKbiYQnxPT4hPRTZNX6zE0z+LrbzzVwlF3/rvs6dtNRh9oE8vabHane5Nh9lAz6OOe5U6OyrBRlAv1Jl7pCySBkukXbp7fVdzFdKiFIr8rqmus8RQ711FdU133V617Uab9j02cJKL2von8Pdeyq9CzzuHOgEoKI1FI12rcuZcjTpxKCLFXsuv4wYCpbUglWMVNr+p0su85ebQXdqYQgvVPRfF6WNYZDv5vpUECQpaqbKQyVeYQ2g+zI2qkPbFtGuofo4ltnX+pavDwKCFJpCqOBx5imIfVxvkNntG1+N6v2PY5N1IBgZlcC5wP3uvtZMdMi/YndTlDHXHqqZNtElllS6CIbLC647NrKgWxNjzvV7zSm2CWEDwLvAbZHToeUWHaRvclxupZepraoT/Z8j7Du+xzq+2zyu3nokfUN3ZKVmccRNSC4+w1mdkbMNEh/6mbOY6yCmmrGczDpJBjGVQxxHkN8P2E8QxgFfbDHzo5j/H1NRewSQiUz2wJsAVhbW4ucmtXQV7/0tsdtciPXeQrNey87UC5WptH0uGG7bFXRmNU9tzNPOYFde/dzzNFHKfOOZPQBwd23AdtgMQ4hcnKkRN0qgqLtyqpvqvbZx8pcU5Ktby+6Ln3OFVRntHITeb+DEOi6UK+k9kYfEGR5ll3UzmbiXW7kvAbJovMJwtP1kNUtRbpe6zpBtGlaYmacu/cdYOv2HdFLbKtOAUF616SKoEze5GpFpYqp1Rf3VaKpKhl0uS55AXrr9h1s3b6j15XQ0iWOPo39NzBGsbudfgx4AXCSmd0JvNndPxAzTatsWUXtqsyq7Om3SUaabrjMm4VzmcEje6y6VT7LSNMQwbRpwFM1zzjE7mX02zGPL8vRtitoesrlkKlnhUw/3f2yz8yk7wyqaNbVvo/XZwab/mzVfsPfVfUzTaoyksMMffO2yayK+qnD4mk0dFsMddFBOqhs3b4jN3MbUlUAGHpVsbIn9SGeyrPHa1tSkDgUEGQwbaskiuqu0/uCRzObdNfLvpZwHKo6paiqaFnH60tRusLxwnloGc1pUUCQaPpo4ISNQQIWmVAICkNOldxlnESbRtShRnAPcW366rmlKqflUkCQwXStkmhSkkiXCvrotTJ0I2fRoL+5ZIAqGUyTAoJMSlVVRWhDyC6m0lVZyaRJSaHKrr2Hr1081AjuvGO0DUTZNHa97lPtTjx1CggyuGU1UkP/A7X6HABWxzFH93NLdklvm8y3z5Hi6RKfLJcCgsxKX1UvRb2V+p6+IcjOTRReX/W6F5WeU91pQorkPYmHwYBN9D2uIl3iyxtDIsNQQJBZ6ePJOGSOffRWWrYmVS15mX/o2nvLnv2H1igoy4yrutV2SXs6TdkuwzIMBQSZpS4lg+DBh9cPZXB5vZX6zKDC0p7pkgHkV1+l/95kgZk86Ub47DQhDz1SPZNqtmpniOA5pYA8dQoIMgt9zd2TnRYbaFWF0jYNXdWpMssLJg89ss5ZayduyOAPevnTedhP+tr0OftpH/uUZhQQZNLqVBE1nWjtgsuu5QiDs9ZO3PD3ZciWDMrmewrzND348PqGEdxVac1rsA1rEJx36dUb/h5GiGelA+cte/Z3qiqS8dC3KIdMuWtfn/34Q7XJrr37D02JUbXfrt1S29i970BhtU7ZsdINtvBoQ/Z5l1592Mpl6fPPjvvILs5TpxTVZxdd6Z8CgoxeWe+aorUOyrapO2V03bR0UVQ9UifAlfXsqXPNuqzHnG5rKJpNVqZHAUFmNQioa5qzdehHGBuWdExn4HkNsemeOUNdx7569qSrxODR6rKHHlnfsFZzNrPPlhTq9ETKW9ei7LNT/g1OmQKCRNP0ib1swFjZHD1NJtUrqg7J/r3O4Kmy94vWRM5rK2iiyTXLk54HKryuqn7KC37K0KdJAUFW9ibOO990UAjy6sbTwWH3vgMcYYsSRd703H0LaQyll6q6+6p+/EVP/9n3yj5XdtzsuhZBaJCO0f4i+RQQZOnq3vRNSgFl2mYmIXPftXc/xxx91KHXVXXv2cbavAwvu55zyNTrpDVdrRMy1aLjlE30V3b9616zvC6s6fPTFBTTooAgh6xKXW7Tp9CDvnGQVvapPBu4sr1/+p7nJ13HnxWeykMa0k/oy5qPKS1UQdVdwnRVS6tjoYAgS9embr/J9l3ltSWEQVvpjK0ogy3LsLucS0hTet/paqN0aSDbFtBXGuruq04JRcZHAUFylT1Fj+XprW066lZFZXvCVI2gDdtkg8FDj6z3Notpkbz6+iB2t9A6bRdN3pPhKCBING3rqZs+cXYJHOljVX2+qOdRXsNvmwwvPYo6BJ1019G6VVNNrkfTUlzRa5kGBQTJVVT0z6uTXvbNX9Zo28QQDdLpvvawqNIJ01Gk1ZlrqGgwWnZ1uPTnsvMuKWOWJhQQZPSK6uzLRhSnt286UrlNJpoNnGHls+zgr7qy6zEE6cFoITBkz7No2u4m10PdP1eTAoKUGmMPkLGko0y6oTetTttMnfUYysYeNOnCKpJm7gXdIUZo8+bNvnPnztjJWHljyYjrpqPpiOi8tQ+q9p2dRiJUZxWtLZx3rFCyyOuhlN1PXtrK1lLIS3OfbQh1jeW3s2rM7CZ331y1nUoI0thYbubY6cirsy8SnvbLMuu8uYSK9lU28rgOrUImeQpLCGZ2DfBqd//WUlNUQiUEGVLTp9eykcHZQJEtEQTZ6q/0TKTprqrphumiPv9VpZRsA3STc+2qSylMuuujhPBB4Doz+xBwmbv/sK/EiXQ1ZIbWtIopveJYVrYraGhs/vwlLy7dLm/cQtfut3X20ef1VPXQ9BQGBHf/pJldDbwJ2GlmHwYOpt6/ouvBzexc4E+BI4H/7O7v7LpPkSp15k5qU6VSNDgs3RW1aIBa6E4aBpEFedNFZ9UdLRxzxPAUOgJIdRvCD4EHgU3A8aQCQldmdiTwXuCXgDuBL5vZZ939a30dQ+ZnyO6Qeb18skEhfbzsRG7Z5SbTo5zD0pRFs6Gmg0hRxt/HOgtVXXizo7PbUJfV6SoMCMnT+xXAZ4FnuvsPej72s4Db3f2byfE+Dvw6oIAgg6iTUaUz9AcfXq8sKaQnsUsvKgMbB46FkkF2/YPsMcM+d+87wFWve9FhGXedkkL6ddFgwlgUFMatrIRwMfASd791oGM/Abgj9fpO4NkDHUtmYhlVD2Hheni0T39eMAnVPNleQeF1eq7/ot5FYfbS9HQUYR9ZRSOU28jrirpr7/5DE/pl1ylos2+VDKanrA3h5wc+dt6s8od1eTKzLcAWgLW1tYGTJHNWJ6NKL0BfZ2bTPOm1E/JkSx3Z0kN6Yfvs0pNAo8xambM0EXMcwp3AaanXTwTuym7k7tuAbbDodrqcpMnYDdEbpqhraPp42cCQzfirJrELbQ+79u7fUJVz9uknHtYYHUoEeZPWpaet6HotwjmFBXdClVffpRAZv5gB4cvAU8zsScC3gZcCvxMxPbIiYmVUIePNW00s/e9jNx3FQ4+sbyitpKuxsj2R6lDmLHVECwjuvm5mfwBcy6Lb6ZUDtleIFGqzqEv4TLanUfb9tLwn/fQCN6GEkG6DSDdQZz/fVy8eVStJEHXqCne/BrgmZhpk9RT114fDxxJUZZJ11zkI+0ln/NnlJfPaDKqqpET6pLmMRBJlJYOqLqt1Gp1Dxh9KFUWN2kVzGmUbufteCU0lA1FAkMlrOutpyNTPu/TqDT2CiuYmSr9XpKxkkF6bOV1SKNtH2UI4Y6OqpvlQQBCpoc969qoMPhtE0t1Us0EtbD9UZqzMfrUoIMhkNZ0iIWSoocfOQT+8Hj+9bdU+606Cl+4dVLZ9mbrrJS+TpqiYHwUEkQaaTnbXVFEQya6CNnTmu8xZUWU8FBBkstpkjnW7mDZZ6jK7UllQ1AOpSXrTXVCLjlNX30EkOztsn/uWOBQQZNTqZGJjq04pWrCmjmwmm522IrtdnbTUXdWtKB1pqh6aN62pLKNWJ+NpmznVWcWrbN/ZQWnZKS+KpsJIB4qiBuayabaz6zAXpS/sZ9fe/Rt6U7VZM7ruOWkltHHSmsoyaXUaLMfaqFm1YE2Zvs4h21OpbLK9Mnmzoqb/Hns6bemXAoI0MpZMtw916vTLzrNoSus8YZvd+w4cmuo6PTah7LN51T1NR1TDo9Nk9PndadqLeVFAkFFqkll3qVZqW8ee3mddYe2Esn1VZfB120vS1yZUGeWdY9OMXBn+vCkgSC1jrJ7pq7oib92DJudV1eCdfpIPsmMSsucSgkf6s2EW1PR0Ftl2iLalna6qqpZkGhQQZNSadCXNU9WfPruWMFRPIJfX7TQvUw7v79p7eC+hhx5ZP7QwTvocshPtZT+TDgZVpYXsOgdlXWebTv+hjH6eFBCkliHriuuO+C3KcMvmIKorvRRml+Uj8xxz9MYn+7LeRSGTTzcCpwenhXM+IllvMNuOkQ1wy54dtUmgUXAZHwUEmbWqQBaqZLK9cMoy0+wUGNnG4aLBb3UmqwvHzBt7AItunXVLMulptbPnXtUTqqgKaExVhtI/BQRpZIiSQdWU0lXvh0yvaa+btPQi9+kn+C5jHEJjbjqNbfaVzvjTJYi8ksyyMuqi6qomc0AVVblJPAoIshLKBpbBIhiEoFAnc6qT8WV79tSpn08HETi8CqnuqOd0usqOWzf4atzBatBIZYmuaRtCmfR8P1Wfy44ADvXyn7/kxbVGMeelLTsgDMrbDLL7KqqqyhshnDdLa1G6qhSdbxD2UzZquq509VmX/Uh9GqksK2fr9h2H9d5JvwcbM5zQINtl0rg+MrBsELllz/5Dg8iKBr8VPam3qesva5TPHqeqjUOmTQFBoqtbPVMmZKqhi+UFl127obG3aY+btnXyeZPTtQkaZQPJ0nMTde0RVXRtmrTfNDlWep+hgVztB+OhgCCTl1dNkx4RHN7LZp7h322nk26TznDstGwQgcPbMerW4Zdl2GXVW9lrk25fkdWhgCCzkO6jH+rs4dF69rZVHG2fXNMBpw95k9V1mZsoL4imSwphJbmi9o82x9S8R+OngCCjVbfapaiaZuv2HbndR7t2UW1zHlCvXj8Er65VQXklg7xBa3lBNH1tjrDxrTchwzkidgJEipRNBlckHTyqulwOJR1wbtmzn937DnTOVC+/8Llc9boXcfbpJ3LspqM4+/QT+fwlL+bMU05oVRIJ1ynsL33d0mkNYzOqjnHBZdfWrmbKfi/hekl8KiHI6OT1uukyeKmqT/3QQaNoKc2sodJVtd9sEM1e/yaN8dmMfehSl6qd+qWAIJ31fXPmlQzaLvASQ9W0EH3tu49qr6JG57rVddlxCWEiv7PWTiz8TPZ4mg5jPBQQpJEu3Q3rfibb172q8bQoYwmqRt8uS9Ouq7GOH7ZtE8BCFdNQmbyCyLAUEKS1oW7O8PnzLr2ag75x+ubQB39ZXUW7qBrVnKfuNeyzeqnLYLbwPZx36dUAG6bnrlJVklImv3wKCFJLWU+VJp+B+jd6GHFcVV1UlbEMMS3C1IJTV1XfXfiu6raXtKWgMawoAcHMXgK8Bfg54FnurgmKJqjrjKBV0tM27Nq7/7CRyOltitSdKjoY4lzqBMa2I4P7SGeTTLYondlpQLqMCld1UDyxSgi7gAuA90c6/qwNcSO1eTJb9tNcUcZSZ4bQuvXlITiFqpEHH17nvEuvnkRJoeh7KDv3ptcyBIJlfdfSrygBwd1vAzCzGIeXni3z5m9aMoBFph3WNi5Ka3b2zbIg1qSePL2Psn1WVXst46m5zj6LqoT6SKeqg+IbfRuCmW0BtgCsra1FTs24LSPzmMJN2qQvfbhGdbu1pvcdZladSsmgaS+s9L/LutDWaU+SaRgsIJjZF4BTct662N0/U3c/7r4N2AaL9RB6Sp5MVN3Mt0lf+rRQJZK3fTZjBRqNpG5SzZZ9Heupuag6qWhwWx/pnMJDx1wNFhDc/Zyh9i35YmceWX2VWLoMumo6O2gT2UVkxqrqd1Hn+tYJkCopTN/oq4ykm7EEh1i6jNwtei/GNY1VMmgazLW2wbTF6nb6G8CfAScDV5vZze4+7orYCRnLDVknA63z3lgaVKeu6fWvu79Vf+iYk1i9jD4NfDrGsVeF+nQPZxWuoTL71aQqIzlkqJu/7dO/MqXh9PnAMMbvRb+ZdhQQZkqZ6bhN5XvpO31TOe9VpYAgUaqXmgQsZR79m+sDg6pKu1FAmDndCOOyqhlWlxlV535txkQBQSbZnVKZRT/mdv3mWvJZFgUEmayq+YnGaFUzrHCeYS6qNjOqrsq1ikkBQQ6Zwg2XzSymGBTmYswZ9RjTNAUKCDI5TWcyHaMppbUP2UkE28z82vaYq3atu1BAkElpMpOpDENVOvOlgCBL1ddsmG1mMpV40k/9oYRX9b11LRkoYDWngCCtxbzR6sxkKsMYqmFcGXd8CggjMfebYYintrleq7lKtyPcsmf/4FOlzP2eGoICggDNbh4VyQX6Lxno9xSfAkJkq3Iz6KlN6izHOcTxpD4FhAnqM1NtE5CUuUuf+v496XfZngJCZMvIXMd0g4whDRLXGH8DY7pHYlJAmJAhG2bb7GvVbx5pp+i31lfJYO7Vr0NSQBiJIUsGukFE8uke2UgBYUKGrF5a1RtAlmfozLeP+yM9LcoqUkCYMTX+yhCW/Xsa8nhD93SaGgWECVLGLlO0rAeULiWDOpPvzZkCwoDG8qOKfXyZh2XXty/zeGGCxHCsVaWAMFFjCTYiTfU5fmZM7Q9zoIAwAPVckDlqm2m2/f2PPZMea7q6UECYGAUbWWXL6qm0qhQQBjD2J5sp0LUbr6Ylg66Z99h+A3N+KFNAmBgFG5krzaMVnwLCgPRjbW7OT1+rZq6Z91zPCyIFBDO7HPhV4BFgN/Bv3P2+GGmZqjn9CGW1NXkImGMmPCaxSgjXA29w93UzexfwBuD1kdIiIzLnp69VNdfvcI6/1SgBwd2vS728EfjNGOmQw83pxy3TUCdjHaIqUb/1w42hDeGVwCeK3jSzLcAWgLW1tWWlSSLTTSpjN8f2LnP3YXZs9gXglJy3Lnb3zyTbXAxsBi7wGgnZvHmz79y5s9+EdjSHHwEc/uM++/QTgemfl8xL1/tt6/Yd7N53gDNPOSH3t95lbfE690ys/MLMbnL3zVXbDVZCcPdzyt43s1cA5wMvrBMMJI65BDyRvi2rDWGZ92CsXkbnsmhEfr67/yBGGrqaW3Fx2Qugi7TRpWQAGyevO3bTUZx5ygmHSgZbt+8Y7H6eSn4Rqw3hPcAm4HozA7jR3S+KlBbJMZUfsEhsbafbDjOsFolxD8bqZfTkGMft0xDFxTFkuioZyBxV3a9DV/+k9x+CwRgfrsbQy0hGaI59rEViCsHgwYfXuWXP/sp7K8Y9qIDQUZ8lA1XPiAyr6p4a+p5L924aIwUEKaWgJNIn3zAyAAAEgklEQVSPtk/8y7wHFRBGQNUzIqtl974DbN2+Y3T3+hGxEyAiskouv/C5lT2MYhlspPIQxjhSWUSkrlgzAtQdqawSgoiIAGpDEBFZmrG3F6qEICIigEoIkzTWpwsRqWes965KCCIiAqiEMCl1RjSr9CAibamEICIigMYhTFJZyUArnolIlsYhiIhIIyohzIzaEEQkSyUEERFpRL2MZkYlAxFpSyUEEREBFBBERCShgCAiIoACgoiIJBQQREQEUEAQEZGEAoKIiAATG6lsZt8B9jT4yEnAdwdKzhjo/KZrzucGOr+xOd3dT67aaFIBoSkz21lnuPZU6fyma87nBjq/qVKVkYiIAAoIIiKSmHtA2BY7AQPT+U3XnM8NdH6TNOs2BBERqW/uJQQREalJAUFERIAVCAhm9nYzu8XMbjaz68zs1Nhp6pOZXW5mX0/O8dNm9rjYaeqLmb3EzG41s4NmNpsufmZ2rpl9w8xuN7N/Hzs9fTKzK83sXjPbFTstQzCz08zsf5jZbclv8w9jp6lPsw8IwOXufra7Px34HPCm2Anq2fXAWe5+NvB/gDdETk+fdgEXADfETkhfzOxI4L3AecDTgN82s6fFTVWvPgicGzsRA1oHXuvuPwc8B/j9OX1/sw8I7n4g9fJYYFat6O5+nbuvJy9vBJ4YMz19cvfb3P0bsdPRs2cBt7v7N939EeDjwK9HTlNv3P0GYH/sdAzF3e92968k/74fuA14QtxU9WclltA0s3cAFwL/APzzyMkZ0iuBT8ROhJR6AnBH6vWdwLMjpUU6MLMzgGcAX4qbkv7MIiCY2ReAU3LeutjdP+PuFwMXm9kbgD8A3rzUBHZUdX7JNhezKM5+ZJlp66rOuc2M5fxtVqXWVWBmxwGfAl6TqYWYtFkEBHc/p+amHwWuZmIBoer8zOwVwPnAC31iA0safHdzcSdwWur1E4G7IqVFWjCzx7AIBh9x96tip6dPs29DMLOnpF7+GvD1WGkZgpmdC7we+DV3/0Hs9EilLwNPMbMnmdnRwEuBz0ZOk9RkZgZ8ALjN3a+InZ6+zX6kspl9CngqcJDF1NkXufu346aqP2Z2O7AJ+F7ypxvd/aKISeqNmf0G8GfAycB9wM3u/qK4qerOzH4FeDdwJHClu78jcpJ6Y2YfA17AYnroe4A3u/sHoiaqR2b2POB/Al9lkacAvNHdr4mXqv7MPiCIiEg9s68yEhGRehQQREQEUEAQEZGEAoKIiAAKCCIiklBAEGkpmfny783sxOT145PXp8dOm0gbCggiLbn7HcD7gHcmf3onsM3d98RLlUh7Gocg0kEyjcFNwJXAq4BnJLOYikzOLOYyEonF3X9oZluBvwJ+WcFApkxVRiLdnQfcDZwVOyEiXSggiHRgZk8HfonF6ll/ZGY/FTlJIq0pIIi0lMx8+T4Wc+LvBS4H/kPcVIm0p4Ag0t6rgL3ufn3y+j8CP2tmz4+YJpHW1MtIREQAlRBERCShgCAiIoACgoiIJBQQREQEUEAQEZGEAoKIiAAKCCIikvj/3yrkc5LSNHAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINEE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.minee import MINEE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100       # batch size of data sample\n",
    "ref_batch_factor = 10  # batch size expansion factor for reference sample\n",
    "lr = 1e-4              # learning rate\n",
    "\n",
    "minee_list = []\n",
    "for i in range(rep):\n",
    "    minee_list.append(MINEE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                            batch_size=batch_size,ref_batch_factor=ref_batch_factor,lr=lr))\n",
    "dXY_list = np.zeros((rep,0))\n",
    "dX_list = np.zeros((rep,0))\n",
    "dY_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Previous results loaded.\n"
     ]
    }
   ],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    dX_list = checkpoint['dX_list']\n",
    "    dY_list = checkpoint['dY_list']\n",
    "    minee_state_list = checkpoint['minee_state_list']\n",
    "    for i in range(rep):\n",
    "        minee_list[i].load_state_dict(minee_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VFX6wPHvm55ASAKhhxBQqkGRJoi6IIIIKro2dO1d1rV314YF1LWXn2DDDqgrsooiFkRQgSAdQQIESIAQWkhIz5zfH/dmMjOZJBNSJhnez/PkyS3n3vPO3DvvnDm3iTEGpZRSgSXI3wEopZSqe5rclVIqAGlyV0qpAKTJXSmlApAmd6WUCkCa3JVSKgBpcg8wIjJfRK71cww9RGS5iOSIyC0+lH9URD60hxNFJFdEgus/0sAlIv8Qke/8HYfyH03uTZCIpIlIvp0EM0XkXRFpXsN1JImIEZGQegjxHmC+MSbaGPNyTRY0xmwzxjQ3xpTWQ1wBydu2NMZ8ZIwZVU/1+b0Boaqnyb3pOssY0xzoBwwE/u3neFx1Btb6OwhX9fQlplSjpcm9iTPGZADfAMme80QkSET+LSJbRWS3iLwvIjH27AX2/wP2L4AhInK0iPwsItkiskdEZlRWr4icLSJrReSA3ZLrZU//ERgOvGqvt7uXZbvY9eSIyDwg3mWesxUqIuNFJMVj2dtFZLY9HC4i/xGRbfYvmDdEJNKeN0xE0kXkXhHZBbxrT79HRHaKyA4Rudau6+garO9O+73cKSJXucQVKSLP2e91togsdFl2sIj8ar9XK0VkWBXvawcR+VxEskRki2u3logMEpEUETlox/d8FdvyShFZ6LKsEZEJIrLRft8fF5GjROQ3e30zRSTMLhsnIl/ZMey3hxPseU8CJ7ts31ft6T1FZJ6I7BORDSJyoUvdY0RknV1vhojcVdnrV3XIGKN/TewPSANOs4c7YbWSH7fH5wPX2sNXA6lAV6A58F/gA3teEmCAEJf1fgI8iPWlHwGcVEn93YFDwEggFKsbJhUI84yhkuV/A54HwoFTgBzgQ8+4gCh7XjeXZZcC4+3hF4HZQEsgGvgfMMmeNwwoAZ6264kERgO7gGPsdX9g13V0DdY30X7NY4A8IM6e/5r9ujsCwcCJdr0dgb12+SD7PdsLtPbyvgQBy4CHgTB7u20GTnd53y6zh5sDg6vYllcCC13Gjf3aWtivvxD4wa4jBlgHXGGXbQWcZ79H0cCnwCyXdbltX6AZsB24yt5u/YA9wDH2/J3AyfZwHNDP35+hI+HP7wHo32FsNCu55wIHgK3A60CkPc/5wbM/vBNclusBFNsfQG8J4X1gKpBQTf0PATNdxoOADGCYZwxelk20k2Qzl2kf4yW52+MfAg/bw92wkn0UIFhfMEe5rGcIsMUeHgYUAREu89/BTtb2+NF2XUf7uL58j/drNzDYfv35wHFeXu+92F+oLtPmliVSj+knANs8pt0PvGsPLwAeA+I9ynjblldSMbkPdRlfBtzrMv4c8GIl26wvsN9l3G37AhcBv3gsMwV4xB7eBtwAtPD3Z+dI+tNumabrHGNMrDGmszFmgjEm30uZDljJv8xWrMTetpJ13oOV5JbYXS5XV1LObb3GGAdWy62jD3F3wEoUhzziqszHwMX28CVYLcg8oDVWkl9md3ccAL61p5fJMsYUeNS93WXcddiX9e01xpS4jOdhtaDjsX7pbPISf2fggrJ12us9CWhfSdkOHmUfoHx7XYP1q2m9iCwVkTO9rKMqmS7D+V7GmwOISJSITLG7mA5ifanESuVnMHUGTvCI+x9AO3v+eVi/XLba3XFDahi3Ogx6kCmw7cD64JUpazVn4iURG2N2AdcBiMhJwPcissAYk+plvX3KRkREsLqHMnyIaScQJyLNXBJ8IlbL0pvvgHgR6YuV5G+3p+/BSkjHGOu4gzee69wJJLiMd3IZ9mV9ldkDFABHASs95m3Harlf58N6tmP9UujmbaYxZiNwsYgEAX8HPhORVlT+3h2uO7F+5Z1gjNllv/fLsb748VLfduBnY8zISuJeCowTkVDgZmAm7u+9qgfacg9snwC3i3UAsznwFDDDbn1mAQ6sPlcAROSCsgNnwH6sD7G3UxJnAmNFZIT9gb0Tqw/31+oCMsZsBVKAx0QkzP4SOauK8iXAZ8CzWH3h8+zpDuBN4AURaWPH31FETq+i+pnAVSLSS0SisPq2y+o5nPW5LvsO8Lx9QDTYPqgZjtWtdJaInG5Pj7APziZ4WdUS4KBYB4Ej7fLJIjLQjudSEWlt13fAXqYUL9uylqKxvugOiEhL4BGP+ZkedX0FdBeRy0Qk1P4baL/PYWKdcx9jjCkGDuJ9n1J1TJN7YHsH66DhAmALVuvyXwB218aTwCL7p/RgrFMqF4tILtbBt1uNMVs8V2qM2QBcCryC1Wo9C+vUzCIf47oEq395H1bieL+a8h8DpwGfenSL3It1IPd3u/vge6wWp1fGmG+Al4Gf7OV+s2cVHs76PNwFrMY64LsP60BukDFmOzAOq3slC6uVezdePnvGOrf/LKw+7i1Y7+1bWAc8wTogvNbePi9hHVguqGRb1saLWAeg9wC/Y3VPuXoJON8+k+ZlY0wOMAoYj/WrbhflB7IBLgPS7Pf0Rqx9R9UzMUYf1qGOTGKdvrkGCPf40lCqydOWuzqiiMi5dldBHFbr8n+a2FUg0uSujjQ3YHWPbMLq+73Jv+EoVT+0W0YppQKQttyVUioA+e089/j4eJOUlOSv6pVSqklatmzZHmNM6+rK+S25JyUlkZKSUn1BpZRSTiJS1RXdTtoto5RSAUiTu1JKBSBN7kopFYA0uSulVADS5K6UUgFIk7tSSgUgTe5KKRWANLkrpVQD+HrVTvYf8vWu2LWnyV2pBpZxIL9BP+T1pdRhWJ2eXev1GGN479c0svOL6yCqxmlndj7//PgPJnz0R4PVqcldqQY2dPKPDJ70g7/DqLVXf0zlrFcXsmL7geoLVyFl634emb2WB75YXUeR1b/H/reWZ+eu97l8UYkDsL7YG4omd9Wk5ReVciCvdq3g8VN/Y/SLC+ooIt8U2h/2xuJAXhEHC2rWcl6zw2q178oucK5j98ECtzIvfv8XL8z7y23aF8vT+WVjlnO8oNh66l52nlX/s3PXM+GjZc75g5/6gZe+3+gc37LnED+t3+01puy8Yowx/LnzIEn3fU1K2j5qc+fb695P4csV1iN1563LJOm+r1mdns27i9J47afy56HvP1SEw2FY8FcWGzNznNOLSx3sc/mVtm1fHmsyav9rxxea3FWTNurFn+k7cZ7btILiUvbmFrJ9Xx77DxU5k0dlft+8j/W7cqos4y/7DhXx4BerKSyp/DWs23GQrvd/zY4D+ZQ6DNe9n8KyrftqVE/fifM49tHvvM5bvm0/uw8WkOmRuMtypoiV3PpOnMegp6xfJNl5xTgchhe/38hLP2x0W+72GSu57O0llcby2k+bmLN6l3N818ECXvi+/Ati+H/mc9W0pQx44ntKHYak+77m6W/XszEzh+Mmfsezczfw2bJ0AM5/4zc+XrLN9zcC+HF9prOFPW9dJrdOXwFYiR7grFcXOsveMWMFe3MLOf7xeXR9YA6Xv7OEkS8s4Pr3U1iUuodbpy+n3+PzcP1+OfOVhTgc9X+rdb/dOEyp2tiVXUCQwPZ9FX/mXjjlN1Z59AWnTR7rdT21bfXXlez8Ykodhr8ycwgLCaJfYhwAj3+1ji+WZ/DR4m2kTR7L6vRsznp1IS+N78u4vh0B+OD3rTgMPDJ7LckdYpi3LpM1Gdn8dv8IADZn5XLWKwv59MYT6d2hRZVxOByGoCChqMRBicNB74fnus1PmzyWfYeKmLJgE6UO69fHH1v3c8MH5S3tPbmFDHjie0KCxDmtpNRBfnEpc9dmVqizuNRaz6GikgpfxF+t2uEcnrcuk3+69FnvyS1k2db9APzf/E0c29F61Ozr8zdxzUldytexcif/OKEzAAcLinlj/iZuH9md0GDvbdurp1lJfO1j5c9Gf+uXzV7L/nd5BmOPbV9h+nfrMvluXflrHfaf+W7zH/pyDU+e28frOuuKJnfVJFXVZ+2Z2KuyP69+D+IVlzpwGEN4SDC/b97L8YmxXssd95h7q7nsy6iotLz75r1f03hk9loAbp2+gnYtIjihayvAagXOW5fJPDuhlLUUjTGc+tzPAIx5+Rc+uvYEVqVnc9Owo3A4DBt359KjXbSzjuveTyG5Y0yF1naZ3zbt5eI3f3ebNmWBe+Ira+GXuLROj37wG6/rA3h3URoAy7cdoOdD7s/ivvnj5W6xebpwym/O4RQ70VeIefNesnIKefrb9SxK3cPO7AJKjWH7vjx6t2/Bgo17aBkVxnMXHsclLq/tmEfKv9ie+PrPSuO/5r2a3912Z3ZB9YVqSZO7qrHv12Uy/6/dPHFOH/7KzKFDbCTNw33blQqKS5n2axrXndyVYJeWHVit6OJSQ+vo8ArLGWN4Z1EaFw3s5OzjrYknvlrHg2N7ISLVF/Yix+6P/mjxNv7WvTXNw0Po1DKK3MISgkWIDAt2K3/RlN84pXtr/rdyB+t35TDnlpMZP/V3rhjS2a3cruwCHv9qndc6v1m9k69X7XSOlyV2Zx1Tf+ehM3vj8NJ9v+tgATkFxXy02L1L4h9vLQZgUJc4zvs/KzE+OKaXc/4P63fzQyX92UCFxO7N+7/6dEdaJs35kyuHJvHLxj1e5y9K9T69Mm8v3OJ1GGDgk9+7jU/52fpCcu3++faRXQQSTe6qxq61W1CPj0tm1AsLGNy1JdOvH+LTsi98/xdTft5MfPNwzu+f4DavrO/cWxfKaz+l8p/v/mLu2l0s2VKz/mSAtxZuYUBSS0YntyMlbR+xUWGc9vzPbmU+XryN5I4tODahYuu6j0t/9ORv1jvjTH5kLtERIax+9HQueONXlqbtZ+YNQ1i8ZR+LXeIs6/5577fyxLcmI5szXynvv3W1O6eAm3w4ba6yLwbPmD2VJXaAJ+dU3io9HDNStvtUbsqCzRVa/a7KvojU4dHkrg5bqf2z21uy3ZmdT9voCII8Wuc5BSUA3PXpSjrFRdrdCtX7z3d/VVqXr278cBmPn5PMQ7PWeJ1fdirewnuHExocxMbMXLbsPcRlgzt7LX+j3c+cU1DChl05LE2zugVcuwrKpO+veGygssQOMOjJpn+qpKrcj1X8OqoreraM8slbv2xmusdZB2VdqkF2V0dJqYO5a3exfV8eQyb9yCs/pmKM4fNl6eQXWQfKXE9Lu2jq717PGpi+ZBtb9x4qr8eHMwuS7vuaktLqTy+sLLG7OunpnzjhqR+49O3FPDRrDavSvZ/H/e3a8p/xp1dzKuU9n6+qtl6l6pLU5hzQ2hgwYIDRx+w1DbtzCpwtybTJY0m672sABiW1ZEma1ZL+99he/LRhN4tS9zJ+YCemL/XtpznA+1cP4pMl2/hmTXmyjI0KZcXDo1iUuodZyzP41D61rTrXndyFN3/ZUn1BpfyssjO4qiMiy4wxA6orV223jIi8A5wJ7DbGJHuZL8BLwBggD7jSGNNw19iqOldS6uDjJds46eh48opK2ZSV67VcWWIH97MJapLYAS5/p+I5zwfyivlieTq3z1hZo3VpYlfK4kuf+zTgVeD9SuafAXSz/04A/s/+r5qYbXvzWL/rIGsysnn5x1Tn9JfG93UOT/qmbg++VaWmiV0pVa7a5G6MWSAiSVUUGQe8b6z+nd9FJFZE2htjdlaxjPKz7fvy+H3zXgYktaRLfDMATnv+Z7fzqssUl5Z33ZWdQqaUatzq4myZjoDr7/B0e1qF5C4i1wPXAyQmJtZB1aomlqbt44I3fqvQL71l0hjmb8jymtjBOrNFKdW01EVy93ZViNejtMaYqcBUsA6o1kHdypa6O4fTnl/A3NtOoUe7aK59L4W/9WjNZYM7s3DjHi59u/ycYc9+6S73z2nocJVS9awukns60MllPAHYUUlZVU/KrrT7atUOVmc04/s/M/n+z0wiQ4O15a3UEaguznOfDVwulsFAtva3173CklLumLGCHS73g/5yRQZJ933NntxC8l1uuPTR4vKrIDWxK3Vk8uVUyE+AYUC8iKQDjwChAMaYN4A5WKdBpmKdCnlVfQV7JPtsWTr/XZ5BbmEJk/7eh9zCEj6wL2W/69OVzN9g3R97V3aB8ypQpdSRy5ezZS6uZr4B/llnESmnZVv30bZFBAlxUTz4hXVlpQGGPv0jBcXlBz/LEjvg88U+SqnApveWaYRKSh1k5hQ6b+7keiXbvHUV74etlFKeNLk3Ijuz83nm2w3ERYXxzqLyM1oaywMllFJNhyb3RmTIpB+9Tvd8jJxSSlVH7wrZQH5N3VPh4cFKKVVftOXeQC55azHtYyKcz7W8cMpvLNmyj8fPSWZQUksqeZyjUkodFk3uDWhndgFb9hxi0+5c50MnnvhqHYUl1d+HXCmlakKTez36cX0mV09L4fs7TnFOG+7xFHRN7Eqp+qCdAfWo7JYAf2zz/iQfpdSRaf5dw+q9Dk3udWDW8gwmfLSs0vn3fKaPWFNKlWsRGVrvdWi3TB24bcYKt/Hv12XSpXUzNmbm+CkipdSRTpN7Pbj2fX02rFKqci2bhdV7HdotU4d25xQw6oWf/R2GUqoWutpPJqtL0eEN347W5F6H7py5kr8yvT9MWqm6duuIbpzXL8HfYbiJCK2flBIb5d5H/fd+HWnXIsKnZZ89/9gK06LDQ7jl1KPrJDaA20/rXuX8pf8+zTl8crf4Oqu3Kprc61BeUWn1hZSqA62jw7l9ZHeeu/A457SOsZE+LfvDnX/jjOR2tar/61tOIjI02C1xpk0ey7J/j2TVo6P474QTCQ7y9pA2ODYhhgnDjnKOv/6PftXW57mm5y/sy+8PjKB72+YAbH5qDO1jvCf7s47rwCUnuD/W871rBnF+/05eyz9/UV+v06ty62ndqpwfERrsHL5jZNVfBHVFk3sdWrZ1v79DUE1UWZIqc/GgiolnwxOjncMzrh/sHD7z2PZ0jI1k4b3DueGUrs7p5/f33qo/qnVzXnBJYCN6tqG5l26DhfcO56NrTyDSJTH9dv+p/PHQSI7pEMOfj4/m7L4d3JZpFh5Ci4hQ+iXGsempMRXWKQKzbz7J7TmcYfbl2af2bMOsfw51Tv/yn0NJahXFu1cNZESvtl5fy3e3/420yWMJquSLBKzE+tS5fZzjaZPH0i8xjsRWUaRNHsvbVwxwK9+3UyzrHx/N5zedyHMXHMfSB0/zXKVXFw3oxGuXWF9UNw8/mrtP78HwHq1Z+9jpbuV6tIv2aX21pQdUD5Mxhj6PfkduoT4Yo6l55eLjOf2YdnT/9zc1Xva4hBhWpmdXWy6pVRRpe/MAKxFfNPV3jmrdjB/uHMZrP6Xy7NwNzrLf3X4K3dtG81dmDqNeWADAnaN68MmS7W7rDA8pT7JdW5d/Gbx6SXnL9/4xvTjruA7ERoWyfmcOn9n3958w7CiaR4Tws33v/4jQYFY+MorU3Tn079ySL1dkcOt066yvf516NG/8vImEuCgS4qK4eFAi7yzawr/H9qJ9jPuvg/CQYFKfPAORypNrmXevGshR8Vbcxs7uHWIiKFvUmPKUf1xCDMd1imX+3cMBGHpUPHeO6s6QST8S3zy82rpq4tSebXjj0v4c06EFh4qsz3NEaDD9O8fRv3McAMcnxrJ82wGmXz+YxJZRnDi54k3+nrZ/xYw9dmyFea6iwhom7Wpy98FTc/5k6oLNfHvbyby7MI2TusXTPCJEE7ufLLrvVIZ6+XCVGdA5jpRKfkXFNw/jrOM6eJ1XZmBSHEvTvC8/+bxjOeOlX3jm/GMrvX5hwd3DSYiLpOsD1oPH4+wzI7wlwBYRIXRva7XkureNZskDI5j/VxbxzcM5p28HZq1wfxzxKxcfX23LL7ljDABtXfqk7xndE4AJw8r7mWMiQ+nfuSUAZx7bgbQ9eVx9UhLREaHcOaqHs9ygLnG8s2iLc72eQqq4MVLa5LEc/cAcShyG4T3aOKcbu+1+6ZDOtLO7U/p0jCG+ufVeDT6qldt6wkKCaB8TyepHRxES5L2+U7q1ZkbKdq/zwOr+aeel60ZEGF1NN9WnNwyhxGGc3Ssnd4tnYFJLTu3ZhhKHqXJZf9Hk7oOpCzYDMPrFXwCq3IFU1bZMGkOX++fUah0dYyP59raT+eKPDKbY28bVZzedyLhXFzpb2BPHHcPo5HYs+GsP5/Xr6Cw3/65h3PjhMtbvcr8e4clz+zhb0ACPn5PMQ7PW0DE2kl7tW7D8oZHENQujdfNwrpq21FkuqVUUzcJDaB8b4dZNUNYgLZvSx06SbaLDee/qQW51t2kRwYUDrC6Zyecdy2VDkmgXE8GubOuOotV9MbkKrcHd6IKDpNJ+49HJ7Un592mH3WJecM9wdmbnu02b8Lejycwu4LLBnYmOCOXrW06iZ7sWBAcJP989jIS4KK/rio6o/OKfx89J5uZTj+bkZ35iVO+2fOfxYJsxfdofVvxgfYG5/HDig2tOOOx1NRRN7spvVj48iuMmfldh+qL7TqVdiwiKShxEhgWTdN/XFcr0bNeC+8e0YEyf9ny7dhf/N3+T2/yyJBAeEsTlQ5KAin3QSfHN+Pa2U1i34yBjXv7FOT3MIyleNrgzZ/ZpT7h9JkhZS3x4zzakTR5LYUkpk+as5/aR3YlxufLwq3+dRHRECB1jIxnVuy23jLCS5yndW/Pb/adW6OLwVNY1AL4fLK0vtekK6RAbSQeP+GOiQnlx/PHO8WM6lP8q6Nzq8E5FDAsJolPLKP54aCTNw0MOq9stkOgBVVWn3rzcOjjVxT5XOMTjQJdr10RMVChpk8fy1b9Ock4LCwmiY2wkwUFCZFgw1TmuUyyuVfxyj9VHe15/q4U+/+5h1a6jd4cWpE0eS2LLiq3FsgN8cc3CKu0rDQ8J5tGzj3FL7GB1j3Ru1YyQ4CCmXj7ArVujusSuDl/LZmGEhTS+1Lby4VGsfHhUg9WnLXdVqZjIULLzi6std3Sb5qTuts7vH9m7LesfH01EaDDFpQ6+X5fJTR/9UeXyyR1j2PzUGEqNqVFXgqe7RnWnk52gzz0+gXOPr9k54FMv78+0RWkktozivxNO5Mc/d9O3U+xhx6P86+bhRzPqGO9n2fhDTFT930/GVeP7elONxgL7TIUyb10+wGs5z7PQyg46hQYHcUaf9vxyz3DuGd2DlH9XfkpZUJAcdmIfZh+oG3p07S4O6dmuBZPPO5agIKFfYhx3nd6j+oUaubIDlEeiu07vwbEJR+6Xs7bcVaViokL5/KYhrMk4SL/EOPokuJ8tkRAXSfr+fIKqOQ2uU8sot7M0Dte3t51MXFTFZDUwqSVpk6s+/exINfe2U9iTqw9YPxJpcndhjOGhL9cwqnc72sVE8OTXf7LvUOB8MG47rRvBIjw37y+36Y+fk0ywCMN7tmbl9mxu/HCZ8wrC/p1bOk+X83TZ4M5M+mZ9vcc97aqBhIcE07Ndi3qvK9C0ah5Oqzo+L1w1DZrcXZQ4DB/+vo0Pf9/m71DqXPuYCG47rTubs3IrJPchXVtydJtou1xkla3gJQ+MIDw0mJjIUPbmFjLpm/WM69uRhK37anSaXk0Mczk/WinlG03uLkpKG+fFCDVx/Sldnefluyo7x9f1eovp1w9mcNdWFcpWpY3LhTGtmoez/vHRhIcEIXJUFUu5Cw8J4s5RDXN/DaWOVJrcXUz65k9/h1BrD4zp5TW5lyk7wHbHyO41TuzeuN4QyVcbnjij1vUqpaqmyR1I359H6u5c3v9tq79DOWyLHxjh9WCjp9ioMP6cOLrebs2qlGocfEruIjIaeAkIBt4yxkz2mJ8IvAfE2mXuM8bU7hrzBjTqhQVN8na9kaHB5BeXclqvtm73EfHG5Z5MPl0cpJRq2qpN7iISDLwGjATSgaUiMtsYs86l2L+BmcaY/xOR3sAcIKke4q0XTTGxA3RqGcl3t/+t2nJHt2nOVUOT6j8gpVSj4UvLfRCQaozZDCAi04FxgGtyN0DZeWoxgPut7BqxL1dk+DuEw/LY2ccwsrf3q+9+vnsYi7fs457PVvHulQMZ3lPPNlHqSONLcu8IuN4GMR3wvCXao8B3IvIvoBng9VJEEbkeuB4gMTHRW5EGU+owCDjvYd3UXHFiUqXzOrdqRudWzTi/X0KVDzFQSgUuX5K7t+zgec7gxcA0Y8xzIjIE+EBEko0xDreFjJkKTAUYMGCAX887POqBOfXyINza6Nq6GZuzDlU6//pTunLFiUnkF/l2H3lN7EoduXw5ZSIdcH3mVwIVu12uAWYCGGN+AyKAhnkKbC1s3lN5Im0Ifz++o9t4r0quwPzwmhPo3CqKW0Z0o2NspPOCI6WUqowvyX0p0E1EuohIGDAemO1RZhswAkBEemEl96y6DLSu7MouoKTUUX3BBuD5IN6n/t6HR8/qzZrHTueVi617Xc+7/RRO6hbPz3cP9/qcS6WU8qbabGGMKRGRm4G5WKc5vmOMWSsiE4EUY8xs4E7gTRG5HavL5krj+kDERuBgQTFvLdjMyz+m+jsUr3q2iyYmMpQrh3YBrCfu1Nfl/EqpwOdTU9A+Z32Ox7SHXYbXAUM9l2tMHpq1hi9XNN6TeFyf+q6UUrV1xFym2Nju7nj5kM4Azie/ez7aTSmlauOIyCjb9+WxJiPbL3WHhwTx7PnHuk17aXxfHhjTC4Bx2vWilKoH4q+u8QEDBpiUlJR6r2dTVi4jnvu53usp0yEmgl/vH8FL32+kf+c4TupWftJQ2YOeXW+pW1zqILegxPnQZaWUqoqILDPGeH8smouAP/3i0rcWN1hdF/RP4Cr7gOitp3XzaZnQ4CBN7EqpOhfwyX1ndkGD1fXsBcdVOX/xAyNq9QBopZTyVcBmmr25heQUFNd7PT3bWRcUvXFp/2rLtm0RQUttpSulGkDAttz7P/F9g9Tz7W2nNEg9SilVEwHZci91NKrrp5RSqsEFZMv9l431e+eDLZPG8H8/b6Jvp9h6rUcppQ5XQCb3yd+sr9e701nWAAAgAElEQVT1iwgThh1dr3UopVRtBGS3zPpdOXW6vpfG962+kFJKNSIBl9zr446P4/p25L8TTqzz9SqlVH0JuOT+f/M31ct6+yXG1ct6lVKqPgRccs84kF9v6w7RJxsppZqIgDugWte3ymkWFuwc/uHOv/FXZm7dVqCUUvUg4JL7xt21O5gaHCRcPKgTeYWlPHluH+cteaH8wdNKKdXYBVxy/2PbgVot//yFxzGub8fqCyqlVCMWcMm9NlxvxauUUk1ZwB1QVUoppcmdPyeO9ncISilV5474bpnIsGC+/OdQ8otL/R2KUkrVmSM+uQMcpzcAU0oFmIBK7tn5vj2c4+tbTmLW8gw6xkbWc0RKKeUfAZXccwtLqi3z/tWDOKZDDMd0iGmAiJRSyj8C6oDq6vTsasuc3C2+ASJRSin/CqjkfufMFdWWEdH7wyilAl9AJfdDRd7PeBnQWe/oqJQ6sviU3EVktIhsEJFUEbmvkjIXisg6EVkrIh/XbZiH5/Rj2hLfPJwZNwzxdyhKKdWgqj2gKiLBwGvASCAdWCois40x61zKdAPuB4YaY/aLSJv6Crgmplw2wDn8zpUDaBER6sdolFKq4fhytswgINUYsxlARKYD44B1LmWuA14zxuwHMMbsrutAa6pFhPtLO7VnWz9FopRSDc+XbpmOwHaX8XR7mqvuQHcRWSQiv4uI36/pf/ni4/0dglJK+Y0vLXdvp5d4PhIjBOgGDAMSgF9EJNkY43b/XRG5HrgeIDExscbB1kR88/B6Xb9SSjVmvrTc04FOLuMJwA4vZb40xhQbY7YAG7CSvRtjzFRjzABjzIDWrVsfbsw+iY3S/nWl1JHLl+S+FOgmIl1EJAwYD8z2KDMLGA4gIvFY3TSb6zLQmkqIi/Jn9Uop5VfVJndjTAlwMzAX+BOYaYxZKyITReRsu9hcYK+IrAN+Au42xuytr6CVUkpVzad7yxhj5gBzPKY97DJsgDvsP6WUUn4WUFeoKqWUsgRkch/Xt4O/Q1BKKb8KyOR+0YBO1RdSSqkAFpDJfVCXlv4OQSml/CpgkvuXKzKcwyHBAfOylFLqsARMFpy1PKP6QkopdYQImOReXOp5RwSllDpyBUxyLyzx/qAOpZQ6EgVMcs8v1uSulFJlAie5V/KIPaWUOhIFTHLflHXI3yEopVSjETDJXSmlVDlN7kopFYA0uSulVADS5K6UUgFIk7tSSgUgTe5KKRWANLkrpVQA0uSulFIBKOCS+y0juvk7BKWU8ruAS+6RocH+DkEppfwu4JL7wKQ4f4eglFJ+F3DJvWf7Fv4OQSml/C7gknvz8BB/h6CUUn4XEMm91KFPYVJKKVcBkdz35hb6OwSllGpUAiK5Zx7U5K6UUq4CIrkX6PNTlVLKjU/JXURGi8gGEUkVkfuqKHe+iBgRGVB3IVZvzuqdDVmdUko1etUmdxEJBl4DzgB6AxeLSG8v5aKBW4DFdR1kdd5dlNbQVSqlVKPmS8t9EJBqjNlsjCkCpgPjvJR7HHgGKKjD+JRSSh0GX5J7R2C7y3i6Pc1JRI4HOhljvqpqRSJyvYikiEhKVlZWjYNVSinlG1+Su3iZ5jyxXESCgBeAO6tbkTFmqjFmgDFmQOvWrX2PUimlVI34ktzTgU4u4wnADpfxaCAZmC8iacBgYHZDH1RVSilVzpfkvhToJiJdRCQMGA/MLptpjMk2xsQbY5KMMUnA78DZxpiUeolYKaVUtapN7saYEuBmYC7wJzDTGLNWRCaKyNn1HaBSSqma8+kuW8aYOcAcj2kPV1J2WO3DUkopVRsBcYVqmTcu7efvEJRSqlFo8snd4XJHyHYxkX6MRKk68Oap8Mkl/o7Cdwe2wfynwXi5M+uhvVBa0vAxKSAAkvvaHQedw9ERfryXe9YGSFtYcXpJobXzlzTym5sV5sCu1TVfLuMPWP1Z3cfjq7SFkLMLFk/x/v43tIJs74nOm52r4Pc33KdlLIMNX1vDefvgy3/CozFwYHvF5WvDGHA43KelLYT9aVUv53DAhm/KX+P0S2D+U7D+K/j82vJkXlIEz3aF/91aef2LXra2HVjLFR1yL1NcAEV5vr8mhwN+fQUKDlZftkxpsbXv1ORLyFEKvzxvbZf8A74v18CafHK/+r2lzuGjWjevm5Xu3VRxx/e0cyX88Hj5+GuDYNpY2LIAUt4tn/7bq9bOv/iNiusA6wP+dBfI3W2NF+ZW3KHXzoLdf1Zc1uGwYnX94BsDf3xg7Xhl6/TFJxfDGye57+TV7bjf3AdvDofPr7HGCw5C5lrf6ktfBms+hx3LYf9W3+MEyE4v3z7TxsJzPeCbe6zhgoNWYtmyAHatqXwdqd9b79Ef71deprQEfpoEq2aWT8tYBm+PshLPgW3W683NglWfWl9ykxNh4fMV12UMrP3CagTs3WRtmyknw7f3lm+nvZvcl/nsalj+oTW83ctdPYyx4tuzseK8/Vthxwr3ZJ23D+bcY70/n10FEz0eSTltLLx0nLW+b+6zkhhYX/xlyXzBs/DJePjiRmu8LCF/fi2s/hSWTLH2VWNvnxUfwuTOsPU397qy1sO8h+DVQeV1P9XBjj3Nem9f7ANPtbfe25xd1ufAUQqLp8Ir/eGZo6z9IG+ftT99MA6++7e1LxRkl9dVXOC+XxcXWEkdrM/nN/fAN3eXz/9rrrUtfnu94vsKVh0/PGYN79tkfUnNmgCZ66A4v2J5hwPWf229h6s+hSfaWdugnonxtZVRxwYMGGBSUmp/tmTSfV87h9Mmj/VtoR3LYcsvMPQWmHkFNGsNY/9jzduzEV61T9G/Yz2U5Fsf3rbHQOFBaGHvgI/GWP8vnw3t+sAzXSrWc+5Uayde+Dy07mkNn/0K9Lvc2tlyd1mJ4mCGVf6yWfDBOdbwPVusdXY5xUpUZU6+C375D1z6OXx4Xvn0E26CobfCS8dCqb3jDJ5gxRzZErqNstbT/0oIjYTlH1hfIgOutj6Iz/e0lulzIZxyF8y5G7b8XL7+0ZOt9Ua2hJiO0KwNvDG0fP55b5cnebBe40m3Q1QrK5Ee2A4FB2DhCzDkZutD5arnmVbr74ZfrIQRGQt5eyHpJNiTCnl7IKw5rJ4Ji16y1nHCjfBicsX3vc0xsNv+knk02/oSbtPLShD7t0BsIsxzOR8guj1cMw82fmfN6zYSti+Bt0eWl7llObx8fPn4dT9aXSiViUmE7G3Q4Xg46Q7ISLHi9tUlM+HjC73Pu+A9OOYc64tp9r+gRUe4bQ1MPcX69XXKPbDgmfLyR42ATT+Uj0e1st5bgIHXWYn09Kfg1f6Vx9NjLJz/NjzZrnzajQutBoE3QaHgKHaf9q8/IHs7vO9x95KB18LStyqv21XX4bD5J9/KgvW65j7gfV7CQGjTG/54zxqP6WTFV5leZ0HiiTD3/vJpx46HVdPdyx1/qZVfSgqtz3h8d9jzl3uZSz6F7qN8fx0uRGSZMaba64gCO7kv/wi+nAChUTDmWesDvnYW/PqyNf+2NeXJ4e9vWQln1OPw3lmVV3jWS3D85RVbPTUxcqJ7clGNy4iH4YeJ/o5CBbL+V1q55DAc2cn90yut/2u/qPX6lVKqXjyaXX0ZL3xN7gHxNOkEyeLj0Cfg0SZ0loFSStWjgEjuC8MrOSKvlFJHqCZ/tkw3Sfd3CEop1eg0+eQ+L/wef4eglFKNTpNP7nUlV4SfIyP8HYZSDepAUBDndmxHWkhA9NAqF5rcbfe3bsXN7dqQERLs71DqxK+REWwIDfV3GKqR+zEqktSwMN6ObeHvUOrEqvAwvm0W5e8wGoUmndwLcqu/9NcBTG4Zx+bQELaHhFDZdWFb7URYKN4ePFU3ioA/w6x60kJCWB0WVqfr3xsUxLqwUEqBG9q14fyE9qREhNOnSyJrw/yf6OdHRnIwqObvb1ZwELOaNwOsR4D5evKuAU5P6MCX9rI18W2zKPp0SSTHY38oBl6OiyG3HveTunJJ+7Y80SqOpRHhbKhk+5e9ip0hIRR7LVF7DuCT6Obki7A9JJg3Ylv4vA0rUwIU2sHniTCreTMM8I8O7bi7TTybQ0Po0yWRR+Nb1rKm2lsaEc6m0Ib/ZdSkk3tRgfdncTuwull2BQeTHhLCRzHRXNW+LWM6deDR+FbOcn+FhvJmjNViKdvZxiV0YGV4mHNnAWtHyhchT4R8Hz/UqaGhvN8imt8jwp0fmrMSOnBhx/ZsDQnhrE4duKRjO6/L/hwZQVpICCkR4TwY39IZx+7gYD6Obs7D8S2dyW5w5wRubms9svDvCe25qGN7drj8xL6qfVsAxndsz6DOCVzQobzOrOAgt6T/UlwMd7SJ57FWcRTUMHkZ4Mp2bfgzLJRiYK3HF9fO4GD+1a41pyQm4Hljh3wRPmrRnJSIcK/rPjUxgYdat2JfUBCvxMVwbJdEXomNwQHc07oVfboksiI8jD5dEt261kqBHaEhPBzfkgIR9gUFkRESXCGJbQkNITvI+igU2/Hc3SYegG0uv362hwRzT5t43oyN4cWWsYD1hTozujkG+Ckqkj5dEtni5YNcDJU2LA4EBfFDlO83vavsi2VhZAQzo5vzVoyVPFdHhDOjRTRXt2/L+R3bOz8TG0ND6dMl0a0rZnFkBOd3bM8bsS3YGRzMnuAgNoeGsDk0hP1BQc5GSXWWhYfzfFwsD8a3dO4D97ZuxVPxLRmb0J4JbdvwWlwsmcEVfyEviozgvRbRrA2ztmXZL88i4NmWsfyrTTwnJiawPyiISzu0ZUBSIgCPxrfkodatWOKy/7wea11B/nl0c05K7EifLonOeYdE+MLl8/1pdDMGdU7givZtWFrJPliVy9q35bk4a39YER5G2Y0O8kQ4GCRc3b4t5yR0YEcD9wo06YuY1m/aTM8PrEvCD4nQzBgOBAVxcucEZ5mhefks8vjgHFtQSF6QkB4SQkFQ5d9vE7P28nDrVhWmr9qyze3BsttCQpgV3YwxuXmcm9De67re2ZnJ1XaidXVm7iGaOxxMbxFNlMPBf3bvYUK7NhXKtSopJcbhYHMlH7L/bd/BWZ2sWyNM25HJlR0q1uUafylwvL3D37rvAJ2Ki7mrbflzbXsWFnF59kGaGUO3omLGdOrAdQeyOfVQPvGlpYxMLH9GeoTD4fY+HlVUxCb7g33NgWzetj9orpKKiknz8loSi4udCXVp2nYGJpU/4fGigznMaBFd6esqi+X1zCyv77WnFqWlTNyzj9vaVv0838ez9jI4v8DtNR9VVMQnOzIZZMd3dk4us6PL7220ess2toSGsCs4mD3BwTxgf1l0LC7hf+nWUypDge0hIYyxt9t32zK4vl0bWjgcdCop4WuXXxynHsrDARwMDuKPCOsL7PttGTzfMpbjCwrZExzMlLiK77OrKIeDvKAgji8oYHlEzY8vfZm+gy7FJQiQHSRsDQ3lHx3aceqhPH6spCvkrJxD/C+64i+nKIeDXoVFLLO/jBOKi0mvpBtxYH4BSys5HnbX3v38p5V1tfjf8vL5uYovyT+2bGN/cDDndmzHQfvL5eXMLG7x2P69Cou4df8BDHCT/Vl8aM8+4kpLuaOKfaVfQYFz24QYQ0klX8IfZ+wiuagIqeeLmJp0cl/+0+f0WnCNWwJoSNGlDnKCm/SPH6WUHzzZ7lTOPr1+bz/QpDNT/La5fBpdR3eCPAya2JVSh2Ndi6p/LdaFJp2dSiWYZ1rV4gZeSinlBx/9NaPe62jSyd0U5vo7BKWUapSadHLP3r/H3yEopVSj1KSTe9+8X/0dglJKNUpNOrkrpZTyrkkn9/QAuVWAUkrVtSad3Pd5ucpNKaVUE0/uk8OreKCvUkodwZp0cl/dKsPfISilVI3dO/Deeq+jSSd3V71b9a4w7YwuZ3Bet/O4tNelPq1jTJcxTDltitu0kKDq7+Y2svNIr+VbRVS8L41ruaQWSW7Tjoo5CoAecT2YPna61+Wu63Mdb456s8p4IoIj+PXiX7nymCvdpndu0ZmbjrsJgEt7Xcr4HuPd5j8x9Ame+9tzbtMmnjjRLWZvOrfoXGU8naIrvz3E1clXc1GPi7iz/50AnNDuBL4c9yUPDX7IrdyozqN4ZMgjvH/G+wztMJQl/1jCuUefC8CpnU5l0smTAOjVshcAx7U+rkJdqy5fxajOo7hrwF38vdvfuXfgvcSGx/L6iNfp16YfAA8PedhrnLHhsW7jpyWexj/7/tM5fmu/W2nXrB3TRk+r9LUOSxjGP3r9g35t+tE8tOorq1dfsdpt/Lxu57mvq9MwAK7ofUWFZWeNm8XKy1ey+JLFVdZR5rURr7H4ksVu++v1x15f7XIfj/mYCcdN8KmO/m3Lf2X/9+z/us17efjLgPX+/HThT16X79yiM3P+PoczupxRZT2LLl7EayNeqzD9tRGvMaHvBCadPIlHhjziNu+HC37g14t/pU1U+T2d7hl4D08MfaLSei7qcRFnH3U2IxJHVFpmbNexlc67oMcFVb2MumGMqfYPGA1sAFKB+7zMvwNYB6wCfgA6V7fO/v37m9pKnpbs/Nufv9/M3TLXfLXpKzNn8xyzLXtbpWWTpyWbwpJCs3jHYnPjvBvNooxF5svUL51lF2xfYJKnJZvTPzvdFJcWm8U7FpsDBQdMQUmB2bR/k8nIyTA7cnaY1P2ppqikyBhjzG87fjPJ05LNtDXTnHUcKjpkjDEmtyjXPPX7UyanMMes37veJE9LNu+tec8YY0x+cb7ZfnC7s+6duTtNSWmJMcaY15e/XiHuz//63BhjTEZOhpm5YaY5VHTInDPrHGd9qftT3V73R+s+MosyFpkZ62eY/OL8Cu9hqaPUuW5XBSUFZu2etcYYY1btXmWSpyWbzEOZbq8vuzDbWX7T/k0meVqy+X7r926vP/NQpvP9OfmTk83TS542i9IXmUUZi8zO3J1udRaXFjuHHQ6H2Zu/1yzesdhsO+i+LV3L7Mnb43wdry9/3WTlZZni0mLjcDjMsl3LzOqs1eaNFW+YSYsneV2HN5sObDJfbPzCTP9zukmelmwm/jrROS+nMMd8u+Vb53hhSaEpLCl0W75sG09aPMlMWzPN7X0qU+ooNTtydph5afPM5399bopKisyOnB3mov9dZFbtXuVcz9bsrW7LZR7KNM+nPG+yC7PNXfPvMvvz91f5Wk7/7PQK+31RSZGZsnKKOVh40Gzav8nttbnuCzmFOcbhcFT7fiVPSzb3LbjPPLLoEZM8Ldnsyt3lnLf5wGaTW5TrLHfLD7dUu76svCxnHH2m9THJ05LNuC/GOd+3/fn7TXFpsZmxfoZJnpZsDhQcMFl5WW7rWLdnnVmwfYFZv3e9OVBwwGs9N8670dy74N4K08vKF5UWmRvm3WDWZK0x+/L3OWNKz0l3Kz/+f+PNV5u+cs4vKClw7hOZhzLNxF8nmg/XfWh25u40e/L2mL35e6t9D6oCpBhf8na1BSAY2AR0BcKAlUBvjzLDgSh7+CZgRnXrrevkXp3ZqbPNVd9eZZKnJZvNBzZXWz4jJ8McLDxYo3hW7l5pSh2lPq3blw+NMVbyyC/ONz9u/dEkT0uukLyNsXZCb4m7PjgcDvPB2g8q/cAEkpW7V5rkacnmmy3f1HjZVbtXOb/4/Sm7MNus37venDPrHLNmz5pqy5c6Sn3ah10VlBQ4GyRVKSop8qmcq8KSQjPui3Hm9x2/12i5+rBk5xK3hpi/+Jrcq70rpIgMAR41xpxuj99vt/gnVVL+eOBVY8zQqtZbF3eFPPnNfhwIK+aEtkN5a/Qb1ZY3xpCVn+X280upqmQXZhMTXvWtdJVqSHV5V8iOwHaX8XR7WmWuAb6pJKjrRSRFRFKysrJ8qLpqIUHWwzoSor0/9MJL/ZrYVY1oYldNlS/J3dsd570290XkUmAA8Ky3+caYqcaYAcaYAa1b1/6Wl4eCrQN8zcJq/hg1pZQKZL482C8dcD3dIQHY4VlIRE4DHgT+ZowprJvwqpYv1ndMVKg+EFcppVz50nJfCnQTkS4iEgaMB2a7FrD72acAZxtjdtd9mFVrFqItd6WUclVtcjfGlAA3A3OBP4GZxpi1IjJRRM62iz0LNAc+FZEVIjK7ktXVC225K6WUO1+6ZTDGzAHmeEx72GX4tDqOq0Y0uSullLuAuEI1KkSTu1JKuQqM5K4td6WUchMQyV0PqCqllLuASO7acldKKXcBkdybhWrLXSmlXAVEco8MifR3CEop1ahocldKqQAUEMk9NCjU3yEopVSjEhDJXcTbvc2UUurIFRDJXSmllDtN7kopFYA0uSulVADS5K6UUgFIk7tSSgWgJpvcHY6qH+ytlFJHsiab3AsO7SfYaIJXSilvmmxyX7X9AEPzC0gsaLIvQSml6k2TzYz7N6VQAsSS7+9QlFKq0Wmyyb1k08+UihCCds0opZSnJpvcjz/4EyVAsOZ2pZSqwKcHZDdGCY4MlkUm+jsMpZRqlJpsy73A6J0glVKqMk02uUdIsb9DUEqpRqtJJne9gEkpparWJJN7vwdmOIcHxvXyYyRKKdU4NcnkviLiBucJkEv3/+nXWJRSqjFqcsndOBwAbAltsif6KKVUvWtyyX3Hju0AjEvo4OdIlFKq8fIpuYvIaBHZICKpInKfl/nhIjLDnr9YRJLqOtAyBzYv44/w8PpavVJKBYRqk7uIBAOvAWcAvYGLRaS3R7FrgP3GmKOBF4Cn6zrQMnsyt3FFh7bO8YcGP1RfVSmlVJPlS8t9EJBqjNlsjCkCpgPjPMqMA96zhz8DRoiI1F2Y5W7Oe9tt/MIeF9ZHNUop1aT5ktw7AttdxtPtaV7LGGNKgGygleeKROR6EUkRkZSsrKzDCvja2DOcw79e/OthrUMppQKdL6eceGuBe15F5EsZjDFTgakAAwYMOKwrkW4d9wy38szhLKqUUkcMX1ru6UAnl/EEYEdlZUQkBIgB9tVFgEoppWrOl+S+FOgmIl1EJAwYD8z2KDMbuMIePh/40Rh9Bp5SSvlLtd0yxpgSEbkZmAsEA+8YY9aKyEQgxRgzG3gb+EBEUrFa7OPrM2illFJV8+kyT2PMHGCOx7SHXYYLgAvqNjSllFKHq8ldoaqUUqp6mtyVUioAaXJXSqkApMldKaUCkPjrjEURyQK2Hubi8cCeOgynrmhcNaNx1VxjjU3jqpnaxNXZGNO6ukJ+S+61ISIpxpgB/o7Dk8ZVMxpXzTXW2DSummmIuLRbRimlApAmd6WUCkBNNblP9XcAldC4akbjqrnGGpvGVTP1HleT7HNXSilVtabacldKKVUFTe5KKRWIjDFN6g8YDWwAUoH76qmOd4DdwBqXaS2BecBG+3+cPV2Al+14VgH9XJa5wi6/EbjCZXp/YLW9zMvY3WPVxNQJ+An4E1gL3NoY4rKXiwCWACvt2B6zp3cBFtv1zADC7Onh9niqPT/JZV3329M3AKfXdrtj3cl0OfBVY4nJXjbNfq9XYN1dtbFsy1isR2Wut/e1If6OC+hhv09lfweB2/wdl73c7Vj7/BrgE6zPQuPYx2pS2N9/WB/UTUBXIAwrmfSuh3pOAfrhntyfKXtzgfuAp+3hMcA39g41GFjs8kHdbP+Ps4fLdr4l9odG7GXP8CGm9mU7KRAN/IX1wHK/xuXyYWpuD4faO+5gYCYw3p7+BnCTPTwBeMMeHg/MsId729s03P6AbLK3+WFvd+AO4GPKk7vfY7LXmwbEe0xrDNvyPeBaezgMK9n7PS6PHLAL6OzvuLAeL7oFiHTZt65sNPtYTROfP//sN3+uy/j9wP31VFcS7sl9A9DeHm4PbLCHpwAXe5YDLgamuEyfYk9rD6x3me5WrgbxfQmMbIRxRQF/ACdgXYEX4rntsJ4NMMQeDrHLief2LCt3uNsd66lhPwCnAl/Zdfg1JpfyaVRM7n7dlkALrGQljSkuj1hGAYsaQ1yUPzu6pb3PfAWc3lj2sabW5+7Lw7rrS1tjzE4A+3+bamKqanq6l+k+E5Ek4HisFnKjiEtEgkVkBVZ31jysFscBYz0w3XN9lT1QvaYxV+dF4B7AYY+3agQxlTHAdyKyTESut6f5e1t2BbKAd0VkuYi8JSLNGkFcrsZjdX/g77iMMRnAf4BtwE6sfWYZjWQfa2rJ3acHcTewymKq6XTfKhNpDnwO3GaMOdhY4jLGlBpj+mK1lgcBvapYX73HJiJnAruNMctcJ/szJg9DjTH9gDOAf4rIKVWUbajYQrC6I//PGHM8cAiru8PfcVmVWY/5PBv4tLqiDRGXiMQB47C6UjoAzbC2Z2XratD3q6kld18e1l1fMkWkPYD9f3c1MVU1PcHL9GqJSChWYv/IGPPfxhKXK2PMAWA+Vl9nrP3AdM/1VfZA9ZrGXJWhwNkikgZMx+qaedHPMTkZY3bY/3cDX2B9Ifp7W6YD6caYxfb4Z1jJ3t9xlTkD+MMYk2mP+zuu04AtxpgsY0wx8F/gRBrJPlZnfdQN8YfVstiM9U1ZdoDhmHqqKwn3PvdncT9484w9PBb3gzdL7Oktsfov4+y/LUBLe95Su2zZwZsxPsQjwPvAix7T/RqXvVxrINYejgR+Ac7EamG5HliaYA//E/cDSzPt4WNwP7C0GeugUq22OzCM8gOqfo8Jq4UX7TL8K9ZZEY1hW/4C9LCHH7Vj8ntc9rLTgasay76PdVxpLdZxJsE6GP2vxrCPGWOaVnK334gxWGeKbAIerKc6PsHqQyvG+va8Bqtv7Aes05t+cNkpBHjNjmc1MMBlPVdjncKU6rFTDsA6dWoT8Cq+nQ52EtZPslWUnxI2xt9x2csdi3W64Sp7+Yft6V2xzkJItXf4cHt6hD2eas/v6rKuB+36N+ByxkJttjvuyd3vMdkxrKT81NEH7emNYVv2BVLsbTkLKwk2hriigL1AjMu0xhDXY1inja4BPsBK0H7fx4wxevsBpZQKRNWGESQAAAA0SURBVE2tz10ppZQPNLkrpVQA0uSulFIBSJO7UkoFIE3uSikVgDS5K6VUANLkrpRSAej/AbwgqkaDNLegAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 10     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            dX_list = np.append(dX_list, np.zeros((rep, 1)), axis=1)\n",
    "            dY_list = np.append(dY_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_list[i].step()\n",
    "                dXY_list[i, -1], dX_list[i, -1], dY_list[i, -1] = minee_list[i].forward()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.plot(dX_list[i, :],label='dX')\n",
    "            plt.plot(dY_list[i, :],label='dY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    minee_state_list = [minee_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'dX_list': dX_list,\n",
    "        'dY_list': dY_list,\n",
    "        'minee_state_list': minee_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current results saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth is 0.4084425498386879 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = mg.ground_truth * d\n",
    "print('Ground truth is {} nats.'.format(mi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "mi_list = (dXY_list-dX_list-dY_list).copy()    # see also the estimate() member function of MINE\n",
    "for i in range(1,dXY_list.shape[1]):\n",
    "    mi_list[:,i] = (1-mi_ma_rate) * mi_list[:,i-1] + mi_ma_rate * mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAELCAYAAADkyZC4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4FFX28PHvIQTCvruxCCgugBAloCAKiiwuwIyiIIqgjrjxw415RUVHYXQcR0dldBxREUQQVFRQEUQURgWRIKgBjawjAZR9DQlZzvtHVZrO3pBUqpI+n+fph1pu3TrpNH1y61bdK6qKMcYYU9oq+R2AMcaYiskSjDHGGE9YgjHGGOMJSzDGGGM8YQnGGGOMJyzBGGOM8YSnCUZE+ohIsoisFZHRRZQbICIqIgnuenMROSQiK93Xf7yM0xhjTOmr7FXFIhIDvAj0BFKAZSIyW1VX5ylXCxgJLM1TxTpVjfcqPmOMMd7ysgXTCVirqutV9TAwHehfQLlxwFNAmoexGGOMKWNeJpjGwKaw9RR3W4iInA00VdWPCji+hYisEJFFInKBh3EaY4zxgGeXyAApYFtoXBoRqQQ8CwwroNxWoJmq7hSRDsAHItJGVfflOoHIcGA4QI0aNTqcccYZpRW7McZEheXLl+9Q1UZe1O1lgkkBmoatNwG2hK3XAtoCC0UE4ARgtoj0U9VEIB1AVZeLyDrgNCAx/ASqOgGYAJCQkKCJibl2mzKyaa/TUG1ap2kxJY0xQSMi//Oqbi8TzDKglYi0ADYDg4DBOTtVdS/QMGddRBYCo1Q1UUQaAbtUNUtEWgKtgPUexmpKYMj7QwBYOGyhv4EYYwLFswSjqpkiMgKYB8QAE1V1lYiMBRJVdXYRh18IjBWRTCALuE1Vd3kVqymZMReO8TsEY0wASUUZrt8ukRljzNETkeWqmuBF3fYkvymx9bvXs363XcE0xuTmZR+MiRI3zboJsD4YY0xulmBMiT3W/TG/QzDGBJAlGFNi3Zp38zsEY0wAWR+MKbHkHckk70j2OwxjTMBYC8aU2K0f3QpYH4wxJjdLMKbEnujxhN8hGGMCyBKMKbEuTbv4HYIxJoCsD8aUWNK2JJK2JfkdhjEmYKwFY0psxJwRgPXBGGNyswRjSuwfPf/hdwjGVEg79qVRr2ZVYiodmf3kmmfm8+CVZxPfoiGZWdmIQEwl52LUsrXbqFujKiNe/SpUft7DlwNwIC2Dq/7xaZnGbwnGlFjHxh39DsGYQPl1xwEqVxIOHc7ijle+ZMrIi4mpJLw0bxV/7h/P+t/3sXnXQS5p1wSAlJ0HeG/pBkZedhYAs77dwL/nHZld/rU7ulE5phIzv1nP3tTD3P/mUnqc1ZgFP24GYHjPM9lz8DBvL16XL5b/e+0rftmytwx+6vxssEtTYit/WwlA/AnxPkdiTMFyvufcuafy+XnzbpofV5u42JjQtjVb9/Lmf9fw2EBnHMhtew+xc38aInDqCXXYdSCduCoxXP30fK7pcgoDOrekciXh86TNvPDJqojimnhHd+b/kMJbX60t4U947D595ArPBru0BGNKrPuk7oD1wZijM37Oj9StXpUbup9WonpUlc9+2MyFrU+kqpsgVJXJC38hppLQ8vjajH1nOQB/uaYDXU4/AYAVG3Yw+s2luep64U9dWZz8G9O+9O8L32tT7+rB+99u4N0l6+l82vE8NqijJZjiJNSqpYkdOuTeeM01cMcdkJoKl12W/6Bhw5zXjh0wYED+/bffDgMHwqZNMGRI/v333Qd9+0JyMtx6a/79Y8bAJZfAypVw99359z/xBHTpAosXw4MP5t//3HMQHw+ffQZ//Wv+/S+/DKefDh9+CM88k3//lCnQtCnMmAEvvZR//7vvQsOGMGmS88przhyoXh3+/W94++38+xcuBGDl3++BL78k/kDNI/uqVYNPPnGWx42DBQtyH9ugAcyc6Sw/8AAsWZJ7f5Mm8OabzvLddzvvYbjTToMJE5zl4cPhl19y74+Pd94/gOuvh5SU3Ps7d4a//c1Zvuoq2Lkz9/4ePeDhh53lSy+FQ4dy77/iChg1ylnu3p18ovyzd3jSZH6tWpdT/zs332dPgduvfoQNu5z3dN6Cf6DAvthq1Mk4xMYaDbj1vJv4U48zeHXBzzQ/sJ1xK2eSUr0+5+x2Jl/8cfJMRr3xTf64KpDzzziBr3/+rcB9k0ZcxLAXvgAgYcd6Ehu2zLX/vn7tqF1J+csHP+ba/kDSh3S/9Lxcnz1ZtMizBONpH4yI9AGex5lw7FVVfbKQcgOAd4CO7nTJiMgDwM04E46NVNV5XsZqjl18TGMITy6mQnvsrD+w+O21wFrm5bkqqsC2uNrcNn01qRnZQE3eqRzHiE5DaH5gB2N/eJ8+Pf4Mu44k7Mktz2dai/zPUr264GcANtZsxJCutwEwc9F43mh5PrPKQXK5+4qz0FdfI2bXTv7Z+lIApnz1H8Z2vo41MbVylR267ksape3n6TbOHyMDzz+Fmy4+I/THzW9xdXgo/ir+tWwKlfv0oUq9y+nVvglnT/0PF//+EwCHK8UQm52FXHMNtL8cUlOZt6CgG3DO8/TnDudZC0ZEYoBfgJ5ACs4Uyteq6uo85WoBHwNVgBHulMmtgbeATsBJwGfAaaqaVdj57BKZf5ZtXgZYZ39ZyFYlMyubKpVjii9ciF0H0qhfMw5VJTNbiY3J/zjcwfQMEtdup1Gdavx39VbeX7qBOtWrcFHbk/jg242hcuNvPp+Rr319zLGUpZpxlTmQlnlUxzx7YxdaHFeLFet3cNbJDUhNzyAtI4vVKbvJzlbGz0ni5VsvpPlxtVidsptd+9NYnPw7tatX4bZerUP1/Lx5Ny2Pr13s7y1blYzM7NClvrLg5YRjXiaYzsCjqtrbXX8AQFX/lqfcczgJZBQwyk0wucqKyDy3rjzXUY6wBOMf64MpO73HfQwcufW0ODv2pbE39TCN61cnrkplJsxfzcxvNvCnHmfwxqJfOJyZDcCs+3uzaefBXLe3eu3Ss5vyyYpNJarj9t6teWneaq7teirfrd9B8pY9APyhU3NuvOh0RIR+T84FnPcs5/0TYO7Dl6Oq9PnrHADeGdWTq5+ez9DupzH4glYliqs8Ka8JZgDQR1X/5K4PAc5V1RFhZc4GxqjqVSKykCMJ5gXgG1V90y33GvCJqr5b2Pkswfgn5yn+tse19TmSiusv05fxzZptofXj61RzrtEn/0af+KZs35fGXZefFdq/ZuveMk0WR+uq81owvGdrvt+4k/83xbnc9a+bzycjK5s2TesDsHnnQZau+Z2X5/+U7/iEUxox7tqOVCrkrrBwO/enUbt6lVBLTVVz3U12IC2DypWEuCqVycjKpnIlKfRus4rIywTjZR9MQb+hUDYTkUrAs8Cwoz02rI7hwHCAZs2aHVOQpuQssRybw5lZPDP7By5qexLnnXY8v+1OZXXKbv63fT9DLzo99OV56HBmruQC8PveQ7y3dAMAkxc6NzjM+e5Xz2L9+5BzuX/KkTuuZo/uE2oZFKTl8bV5afgFwJFWV7+OJzN7mdNJP7ync/moffMGhbbGGjeowR/rt+CskxvQqHYcNeNiSc/IokZc7FHF3qBWXK71vMmjZlh9BV0uNMfOywSTAjQNW28CbAlbrwW0BRa6v/ATgNki0i+CYwFQ1QnABHBaMKUZvInc4k2Lgege9LKg5yz2HEznoWnf0r3NSVSvWpnLO5wc2vfx8v8xfo7T8lu4Kt9Hm+lf539grqw8MbgTSZt2sS/1MB8t/5XBXU8lvnnD0P7j6lSjamwMD1x5NpO+SObV252HAAsTnkBuuviMiFodOUSEVifWCa0XdR4TPF5eIquM08nfA9iM08k/WFULfAIpzyWyNsA0jnTyLwBaWSd/MFXkPphdB9J4Z8l6bu3Zmqxs5bLHnev1s0f3IXnLHv78xjcM7noq09wH5abe1YPPkzazdM02Kgn88L9dpR5TzbhYDqRlHPVxjw1M4Lg61Vi5YQeHM7PpfPrxnFC3OvNWbmLVpt30jm9Ku5PrF/klnpXtfF+ED11iyrdy2QcDICKXAc/h3KY8UVUfF5GxQKKqzs5TdiFugnHXHwJuAjKBu1X1k6LOZQnGPzmzWZ7e8HSfIym5BT+kkLRpN1d3bsmXP21l4ufOz9a2WX2Sfi39ZHE0GtWO44Erz6Z1k3r8uuMAJzdybnVNTc9k94F0kjbt4oNvN9KsYU3aNK3HmU3q0erEOvy6fT93vvoVs0f3iaq+BROZcptgypIlGLP/UAY14iqHLsEkb9nDyY1qhYb/yFZlzda9jHzta6be1YOGteP4efMenvvoB/45rAt/fKrsH7XK25cxa3Qf+oetP39TF+rXjGPb3kO0bVa/zOMzFZ8lmAhYgvHPoo2LAOjWvJvn58rIyiYrW4mLjWHdb/t48v0VPHtjl3yjxHY6tRHfrt0OwF+v7ciYt5Z5HltR6lSvQtczT6Dl8bWpV6MqY99ZzuzRfQp83uFAWgaHDmfSqHY1HyI10cYSTAQswfinLPpgdh9IZ9Czn3lWf6TaNqvPwwPO4U8vLeJgWiazRvcOPTz3xMzvWLR6K3+//lzaNW/Aph0HaNawJhklfDDSGC9ZgomAJRj/rN+9HoCW9VoWUzIymVnZvLrgZxJOacRpJ9UhNqYSf/h72V2+evbGLlQS4a6JzhPqc8dcxncbdvDg1G+ZO+ayQvsxsrKzOZieSe1qVcosVmNKyhJMBCzBlG+l0UL5f/3b83nSFhLXbS9w/yNXdwiNqgtHbp/dsusg+w5lcEbjurnK7009jKpSt0bVEsVlTJBZgomAJRj/fLbeSQyXtLykyHKHM7MY+MxnzLjvEtIzshkyfgGDL2jFa+6ghpF4654e/PXd71i1aTdPDTmPhau2MOe7X5ky8mKOq+P0WRxMz+Dj5b/yx3NbMGXhL/y64wCPDvTk/48x5Z4lmAhYgvFPJH0whzOz6Pu3wp/8jlSkY3AZYyJTXoeKMVFiyh+n5NuWlpHFwqTN9IpvyuadB/nTS4uOqs4/dGrOkG6nsWHbft5dsp7MrGweH9yptEI2xpQBa8GYUpeannlMz5Tc168dvdo3Lb6gMabUWAvGBNrctc6lrz6n9mHZ2m3FPnPy+OBOPDTtW3q2b8Kofu3LIkRjjA8swZgSu3/eY2zedZAOFDpUHAD39m1Hh5aNaFg7jo8evNRGrjWmgrMEY0pkdcpu6u28k3qQb5KF6y5oxZXntaBmXGy+OTgsuRhT8VmCMUftjYW/cN5px/F/7lS5VaVevjJ5h0GxQRaNiT6WYMxR2bEvjalfrmHql2tC27apMxHVcXIugF3+MsYAlmBMhF7+dHVoBsW8dlT5iObH1WLGwIeoVEksuRhjAEswpgjL123nxbmreO2OboUml/bNGzD1Kmck47rVbUgVY8wRlmBMgXbuT+PBad8C0Oevc/LtH3j+KbQ8vjbd25xU1qEZY8oJT69liEgfEUkWkbUiMrqA/beJyI8islJEvhKR1u725iJyyN2+UkT+42WcJr/Bzy0odF+7k+tz08VnhJLLez+9x3s/vVdWoRljygnPWjAiEgO8CPQEUoBlIjJbVVeHFZumqv9xy/cD/gn0cfetU9V4r+IzBdu5P63Q5FLYOGDjl44H4Mozr/QsLmNM+ePlJbJOwFpVXQ8gItOB/kAowajqvrDyNYCKMW5NOTb0X1/kWq9auRIXndWYkZedVegxswbN8josY0w55GWCaQxsCltPAc7NW0hE7gTuBaoAF4ftaiEiK4B9wBhV/bKAY4cDwwGaNWtWepFHKVUlIys7tD6qX3t6tm9S7HF14up4GZYxppzysg+moCfr8rVQVPVFVT0FuB8Y427eCjRT1bNxks80EaldwLETVDVBVRMaNWpUiqFHn7TDmbk68wedf0pEyQVgRtIMZiTN8Co0Y0w55WULJgUIHxq3CbCliPLTgZcAVDUdSHeXl4vIOuA0wIZLLmW/7UnNd1kM4MaLz4i4jpcSXwJgYNuBpRaXMab88zLBLANaiUgLYDMwCBgcXkBEWqlqziPhlwNr3O2NgF2qmiUiLYFWwHoPY41Ka7buZcSrX+Xbfm/fdkdVz5zr8t/GbIwxniUYVc0UkRHAPCAGmKiqq0RkLJCoqrOBESJyCZAB7AaGuodfCIwVkUwgC7hNVXd5FWs0SsvIKjC5HMuMkdVjq5dGSMaYCsYmHItC7y5Zzyuf/ZRr23392nHK8bU55YSj77B/84c3Abi+3fWlEp8xpuzYhGOm1GSr5ksuJZ3n/tXvXgUswRhjcrMEE0W27T3EkPGf59r28IBzSlzv/CHzS1yHMabisQQTBVSVVZt2c9/kJbm2X96hGV3PPLHE9cfGxJa4DmNMxWMJJgoUNFhlSS+LhZu0chIAw+KHlVqdxpjyzxJMBaaqBSaXEZe2KdXzWIIxxhTEEkwFtfjn33jsneX5tn/84KVULuUJwRYOW1iq9RljKgZLMBXQ6pTd+ZLL368/l/bNGyBS0Ag+xhhT+izBVDDDXviCrbtTc20b1a+9p8nlleWvAHBLh1s8qd8YUz5ZgqlAsrI1X3L58IE+VKkc4+l5Z6xyBrq0BGOMCWcJppz7dft+Rr3xDXtTD+faHt+8AX8fcl6ZxPDZDZ+VyXmMMeWLJZhy7pb//DfftsEXnMrQ7qf7EI0xxhxhCaacWrVpF/dOWlLgvrJOLv9e9m8A7uh4R5me1xgTbF5OOGY8cjgzq9DkMu3uHmUcDXz4y4d8+MuHZX5eY0ywWQumHOr7t7n5tpXmk/lH65PrPvHt3MaY4LIEU05s3Z1KtSoxDPyndagbY8oHTxOMiPQBnseZcOxVVX0yz/7bgDtxJhU7AAxX1dXuvgeAm919I1V1npexBlnvcR8XuP3lWy9k36HDtDu5QRlHlNvz3zwPwF3n3eVrHMaYYPGsD0ZEYoAXgUuB1sC1ItI6T7FpqnqWqsYDTwH/dI9tjTPFchugD/Bvt76oszpld4HbZ9x7Cc2Pq+V7cgFYsGEBCzYs8DsMY0zAeNmC6QSsVdX1ACIyHegPrM4poKr7wsrXAHKm1+wPTFfVdGCDiKx16yu4Z7sCu+f1xfm2XdT2JOrWqOpDNAWbfe1sv0MwxgSQlwmmMbApbD0FODdvIRG5E7gXqAJcHHbsN3mObexNmMGU9OuufPO33NrzTH7evIf7/xDvU1TGGBM5LxNMQQNfab4Nqi8CL4rIYGAMMDTSY0VkODAcoFmzZiUKNkj++eH3zFuZkmubn3eJFefpxU8DMKrLKJ8jMcYEiZcJJgVoGrbeBNhSRPnpwEtHc6yqTgAmACQkJORLQOVV3uQSdEtSou7KpTEmAl4mmGVAKxFpAWzG6bQfHF5ARFqp6hp39XIgZ3k2ME1E/gmcBLQCvvUw1sAY+3Zivm1/uaaDD5FEbuY1M/0OwRgTQJ4lGFXNFJERwDyc25QnquoqERkLJKrqbGCEiFwCZAC7cS6P4ZZ7G+eGgEzgTlXN8irWIPk6+ffQcpAvixljTHFEtWJcWUpISNDExPx//ZcXGVnZXPHEkSfiP7i/N9WqlI/nYJ/8ynm8aXTX0T5HYow5WiKyXFUTvKi7fHyDRYHw5AKUm+QCsPK3lX6HYIwJoPLzLVaBvTg3Kdf6zD/38imSYzN9wHS/QzDGBJAlGJ9d88z8XJOF3dmnDTXjYn2MyBhjSoclGJ+FJ5eXb72Q5sfV8jGaYzNu0TgAHu72sM+RGGOCpNixyETkNBFZICJJ7no7ERnjfWgV35Dxn4eWxw5KKJfJBSB5ZzLJO5P9DsMYEzCRtGBeAf4MvAygqj+IyDTgr14GFg227T0UWj631fE+RlIyb175pt8hGGMCKJLRlKurat6HHDO9CCaahA/Bf07Lhj5GYowx3ogkwewQkVNwxwITkQHAVk+jquCWr9+ea33MgHN8iqR0PPLFIzzyxSN+h2GMCZhILpHdiTPe1xkishnYAFznaVQV3INTjzQI5465DJGCxvYsPzbt21R8IWNM1IkkwaiqXiIiNYBKqrrfHV/MlNB1F7Qq98kF4PX+r/sdgjEmgCK5RDYTQFUPqup+d9u73oVUse0/lBFavqH7aT5GYowx3iq0BSMiZ+BMWVxHRK4M21UbiPM6sIooKzubAU9/CkCDWsGZkbKkHvjsAQD+dsnffI7EGBMkRV0iOx24AqgL9A3bvh+4xcugKqpftx8ILf+5f8WZlXLnoZ1+h2CMCaBCE4yqzgJmiUhnVbUZpUrBbRO+DC2f3aLi3Jo8oe8Ev0MwxgRQJJ38K0TkTpzLZaFLY6p6k2dRVUC/7U4NLT89tLOPkRhjTNmIpJN/CnAC0BtYhDN98f4ij3CJSB8RSRaRtSKSb7IQEblXRFaLyA/ucDQnh+3LEpGV7mt2ZD9OMO0/lMHQF74IrZ/VrL6P0ZS+UZ+OYtSno/wOwxgTMJEkmFNV9WHgoKpOxpna+KziDhKRGOBF4FKgNXCtiLTOU2wFkKCq7XDuTHsqbN8hVY13X/0iiDOwcjr2wZlIrKI5lHGIQxmHii9ojIkqkVwiy7mvdo+ItAV+A5pHcFwnYK2qrgcQkelAf5xpkAFQ1S/Cyn8DXB9BveXK73tSc62Xp4nEIvXi5S/6HYIxJoAiacFMEJF6wMPAbJwE8VTRhwDQGAh/xDvF3VaYm4HwaR3jRCRRRL4RkT9EcL5AGj11aWj5kzGX+RiJMcaUrWL/nFbVV93FRUDLo6i7oEfUtcCCItcDCUC3sM3NVHWLiLQEPheRH1V1XZ7jhgPDAZo1a3YUoZWNhUlb2LLLacFMu7sHlSrAU/sFuXvu3QA81+c5nyMxxgRJsQlGROoCN+BcFguVV9WRxRyaAjQNW28CbCmg/kuAh4BuqpoeVv8W99/1IrIQOBvIlWBUdQLOOGkkJCQUmLz88sInSXyY+L/QeoNa9myqMSa6RNIhMAenf+RHIPso6l4GtHLHLdsMDAIGhxcQkbNx5pnpo6rbwrbXA1JVNV1EGgLnE9llucAITy6N69fwMRLvWcvFGFOQSBJMnKree7QVq2qmiIwA5gExwERVXSUiY4FEVZ0N/AOoCbzjDvr4q3vH2JnAyyKSjdNP9KSqri7wROXAxDu7+x2CMcaUOVEt+sqSiNwDHAA+AsIvYe3yNrSjk5CQoImJiX6HAcDlj88hM9t5X+c9fLnP0Xjvzo/vBOxuMmPKIxFZrqoJXtQdSQvmME5L4yGOdNIrR9fhHzVWbdoVSi7RolpsNb9DMMYEUCQJ5l6chy13eB1MRXDvpCPDtn304KU+RlJ2nu71tN8hGGMCKJLnYFYBqcWWMry9ONdNbsTGRPL2GmNMxRRJCyYLWCkiX5C7D6a425SjzmsLfg4tz42ihyqHfzgcsFGVjTG5RZJgPnBfpghZYf0uH9zfu0JMhRypBtUa+B2CMSaAInmSf3JZBFLeXf74nNByRRxvrCg2k6UxpiBFTZn8tqpeIyI/UsAQL+4IyAZY8ENK6A0aN6ijr7EYY0xQFPWn9l3uv1eURSDl2VOzvg8td2p1nI+R+OPGWTcC8Hr/132OxBgTJIXe5qSqW93FO1T1f+Ev4I6yCS/4ssMeVH18cCcfI/FP09pNaVq7afEFjTFRJZLOgp7A/Xm2XVrAtqg0/uMfAahXoyoJpzTyORp/jL1orN8hGGMCqKg+mNtxWiqniMgPYbtqAV97HVh58ckKZ8qbkxvV9DkSY4wJlqJaMNNwJgD7GzA6bPv+oI1D5pfe4z4OLT9xXXReHgO4/j1nItI3r3zT50iMMUFSaIJR1b3AXhEZA/zmDp3fHWgnIm+o6p6yCjKIwpNL9zYnEVMpep/aP73B6X6HYIwJoEj6YGYCCSJyKvAazrTJ04DoeVQ9j7krfs21/sCVZ/sUSTA83O1hv0MwxgRQJH92Z6tqJnAl8Jyq3gOc6G1YwfbsRz+Glj+4v7ePkRhjTHBFkmAyRORanGmTP3K3xUZSuYj0EZFkEVkrIqML2H+viKwWkR9EZIGInBy2b6iIrHFfQyM5X1mb9/DlUffUfkEGvTuIQe8O8jsMY0zARPLteCNwG/C4qm5wp0AutjdXRGKAF3Fuc04BlonI7DwzU64AElQ11b1r7SlgoIjUB/4CJOCMIrDcPXb30fxwXvh4+f+KLxRl4k+I9zsEY0wARTIW2WoRuR9o5q5vAJ6MoO5OwFpVXQ8gItOB/kAowajqF2HlvwGud5d7A/Nz7lYTkflAH+CtCM7rqfFzkvwOIXBGd83XODXGmOIvkYlIX2AlMNddjxeR2RHU3RjYFLae4m4rzM04t0Ufy7Fl7uMomUzMGGOOVSR9MI/itEb2AKjqSqBFBMcVNF59gXMJi8j1OJfD/nE0x4rIcBFJFJHE7du3RxBSySxJ/j20XNkmEwu56u2ruOrtq/wOwxgTMJF8S2a6z8SEi2TS+RQgfICqJsCWvIVE5BLgIaCfqqYfzbGqOkFVE1Q1oVEj74dp+WDZBgB6tW/i+bnKk85NOtO5SWe/wzDGBEwknfxJIjIYiBGRVsBIYHEExy0DWrk3BWwGBgGDwwuIyNnAy0AfVd0Wtmse8ISI1HPXewEPRHBOz6RnZLFyw04AbrnkTD9DCZxRXUb5HYIxJoAiacH8H9AGZ7rkacBe4O7iDnKfnRmBkyx+At5W1VUiMlZE+rnF/gHUBN4RkZU5fTtu5/44nCS1DBjr9/A0/Z6cG1quXb2Kj5EYY0z5EMldZKk4l7AeOtrKVXUOMCfPtkfCli8p4tiJwMSjPacpe/3ecv5emH1tJPd+GGOihT0lGIFla49cvXvl9m4+RhJMPVr08DsEY0wAWYKJwJi3loWWmzW0Yfnzuuu8u4ovZIyJOnavbTF+3nxk8IAZ9xZ6Rc8YY0weRU049i+KuB1ZVUd6ElHA3DXRuWGuTdN61K1R1edogunSqc5Dp59c90kxJY0x0aSoS2SJZRZFOTDu2o5+hxBYfU/r63cIxpgAKmrCscllGUgQpaZnAtDn7KbUqBrRANJR6Y6Od/gdgjEmgIq6RFbkPaeq2q+o/RXBH5+aB8CK9Tt8jsQYY8qfoi4PflfrAAAfqElEQVSRdcYZcPItYCkFjw9WYR1MzwgtX3dhKx8jCb5L3nBufvjshs98jsQYEyRFJZgTcOZyuRZniJePgbdUdVVZBOa3K5/6NLTcO75pESXNwDYD/Q7BGBNARfXBZOEM0T9XRKriJJqFIjJWVf9VVgH64XBmVmh5+j12a3Jxbulwi98hGGMCqMgHLd3EcjlOcmkOjAfe8z4sf/X925Fxx+rVtFuTjTHmWBTVyT8ZaIszCdhjqhp1UzmOuLSt3yGUC90ndQdg4bCFvsZhjAmWolowQ4CDwGnASJFQH78Aqqq1PY7NF6rOs6U14yrTN+Fkn6MpH4bFD/M7BGNMABXVBxOVw8j8ZYbzfOmBtEyfIyk/LMEYYwoSlUmkKEvXOCMn33jR6T5HUn5kZGWQkZVRfEFjTFTxNMGISB8RSRaRtSIyuoD9F4rIdyKSKSID8uzLcichW1ncQ5+lZema30PLg7qeWhanrBB6TulJzyk9/Q7DGBMwng3XLyIxwIs4z9KkAMtEZLaqrg4r9iswDChozt1DqhrvVXwFeWS6Db92LP50zp/8DsEYE0BezgfTCVirqusBRGQ60B8IJRhV3ejuy/YwjojsOZgeWv5kzGU+RlL+XN/uer9DMMYEkJeXyBrjDDWTI8XdFqk4EUkUkW9E5A+lG1p+A/95ZJiTShJVo+KUWGpGKqkZqX6HYYwJGC9bMAV9Sxc6v0wBmqnqFhFpCXwuIj+q6rpcJxAZDgwHaNas2TEHejDtSAf15P+76JjriVaXTXVafPYcjDEmnJcJJgUIH8SrCbAl0oNVdYv773oRWQicDazLU2YCMAEgISHhaJJXLv+Y9X1o+YS61Y+1mqh1e8LtfodgjAkgLxPMMqCViLQANgODcAbNLJaI1ANSVTVdRBoC5wNPeRXokl+cu8eev+l8r05RoQ1sa4NdGmPy86wPRlUzgRHAPOAn4G1VXSUiY0WkH4CIdBSRFOBq4GURyRmp+UwgUUS+B74Ansxz91mpWZJ85NbkMxrX9eIUFd7etL3sTdvrdxjGmIDxsgWDqs4B5uTZ9kjY8jKcS2d5j1sMnOVlbDkefdu5Nblts/plcboKqf/0/oD1wRhjcvM0wQTd9n2HQsuPXpPgYyTl28hzR/odgjEmgKI6wVz//OcAHF+3GrWqxfocTfl15ZlX+h2CMSaAonYsspxRkwEes9ZLiexI3cGO1B1+h2GMCZiobcHsTT0cWm5xfIWceaDMDHjbGUbO+mCMMeGiNsHkPLl/R+/WPkdS/t3X+T6/QzDGBFDUJpgcJzeq5XcI5V7f0/v6HYIxJoCisg9mw+/7QsvxLRr6GEnF8NuB3/jtwG9+h2GMCZiobMHcNuFLwAa1LC2D3h0EWB+MMSa3qEswuw6khZbnPHSpj5FUHKO75ptLzhhjoi/BrNywM7Qs1oIpFX1O7eN3CMaYAIq6Ppi0jCwA3rBh+UvNpr2b2LR3U/EFjTFRJepaMM9//CMADWtX8zmSimPI+0MA64MxxuQWdQkmR0wluzxWWsZcOMbvEIwxARRVCSYr+5jnJDNFuKTlJX6HYIwJoKjqg9l9IB2AwRec6nMkFcv63etZv3u932EYYwImqlow29zh+W1isdJ106ybAOuDMcbk5mkLRkT6iEiyiKwVkXwPS4jIhSLynYhkisiAPPuGisga9zW0NOLZvtdJMMdZB3+peqz7YzzW/TG/wzDGBIxnLRgRiQFeBHoCKcAyEZmdZ+rjX4FhwKg8x9YH/gIkAAosd4/dXZKYclowx9WxBFOaujXv5ncIxpgA8rIF0wlYq6rrVfUwMB3oH15AVTeq6g9Adp5jewPzVXWXm1TmAyV+mm/S58kA1IizycVKU/KOZJJ3JPsdhjEmYLzsg2kMhD99lwKcW4JjG+ctJCLDgeEAzZo1K7bSTLuLzBO3fnQrYH0wxpjcvEwwBT1oEuk3fETHquoEYAJAQkJCsXW3aVrPhofxwBM9nvA7BGNMAHmZYFKApmHrTYAtR3Fs9zzHLixpQKnpmZxQt3pJqzF5dGnaxe8QjDEB5GUfzDKglYi0EJEqwCBgdoTHzgN6iUg9EakH9HK3lcj+tAxqVrP+l9KWtC2JpG1JfodhjAkYz1owqpopIiNwEkMMMFFVV4nIWCBRVWeLSEfgfaAe0FdEHlPVNqq6S0TG4SQpgLGququkMe0/lEEt6+AvdSPmjACsD8YYk5unD1qq6hxgTp5tj4QtL8O5/FXQsROBiaUVS2ZWNukZWXYHmQf+0fMffodgjAmgqHmSPzU9E4AaVaPmRy4zHRt39DsEY0wARc1YZKEEE2cJprSt/G0lK39b6XcYxpiAiZpv24NugqleJWp+5DJz99y7AeuDMcbkFjXftqnpGYA9xe+F5/o853cIxpgAipoEE2rBWB9MqYs/Id7vEIwxARQ137ZpGVkAxMXG+BxJxbNss3M3uXX2mxwZGRmkpKSQlpbmdyjGFRcXR5MmTYiNLburOFGTYA5nOgmmqiWYUvfn+X8GrA/GHJGSkkKtWrVo3ry5Dc8UAKrKzp07SUlJoUWLFmV23qhJMOluC6ZqZUswpe2Fy17wOwQTMGlpaZZcAkREaNCgAdu3by/T80ZRgnFmBLAWTOlre1xbv0MwAWTJJVj8+H1EzXMwoRZMbNT8yGVm8abFLN602O8wjAmUhQsXcsUVV+TbvnLlSubMmVPAEcV74okjI5dv3LiRtm2D/cdd1HzbpmdmUbmSEFMpan7kMvPgggd5cMGDfodhzFHLzMws83MWlWCKiyc8wZQHUXSJLIsqdnnMEy9f8bLfIRiTz7hx45g6dSpNmzalYcOGdOjQgVGjRtG9e3e6dOnC119/Tb9+/RgwYAA33XQT27dvp1GjRrz++us0a9aMYcOGccUVVzBgwAAAatasyYEDB1i4cCGPPvooDRs2JCkpiQ4dOvDmm28iIsydO5e7776bhg0bcs455+SL6fDhwzzyyCMcOnSIr776igceeICffvqJLVu2sHHjRho2bEivXr1ITEzkhRecvs0rrriCUaNGMXfuXA4dOkR8fDxt2rTh8ccfJysri1tuuYXFixfTuHFjZs2aRbVqwZkSPmoSzOHMbLtF2SOnNzzd7xBM0HXvnn/bNdfAHXdAaipcdln+/cOGOa8dO8D9kg9ZuLDI0yUmJjJz5kxWrFhBZmYm55xzDh06dAjt37NnD4sWLQKgb9++3HDDDQwdOpSJEycycuRIPvjggyLrX7FiBatWreKkk07i/PPP5+uvvyYhIYFbbrmFzz//nFNPPZWBAwfmO65KlSqMHTs2VwJ59NFHWb58OV999RXVqlVj0qRJBZ7zySef5IUXXmDlSmdYpo0bN7JmzRreeustXnnlFa655hpmzpzJ9ddfX2TsZSlqrhelZ2RRpXLU/LhlatHGRSzauMjvMIwJ+eqrr+jfvz/VqlWjVq1a9O3bN9f+8C//JUuWMHjwYACGDBnCV199VWz9nTp1okmTJlSqVIn4+Hg2btzIzz//TIsWLWjVqhUiclRf9P369TumlkeLFi2Ij3cedO7QoQMbN2486jq8FDUtmPSMLLuDzCN/WfgXwJ6DMUUoqsVRvXrR+xs2LLbFkpdq0TOo16hRo9B9OXdbVa5cmezs7FB9hw8fDpWpWrVqaDkmJibUd3Ksd2qFxxN+XqDIh1XzxnHo0KFjOr9XPP2TXkT6iEiyiKwVkdEF7K8qIjPc/UtFpLm7vbmIHBKRle7rPyWNJT0zy56B8cjE/hOZ2L/Upu4xpsS6du3Khx9+SFpaGgcOHODjjz8utGyXLl2YPn06AFOnTqVr164ANG/enOXLlwMwa9YsMjIyijznGWecwYYNG1i3bh0Ab731VoHlatWqxf79+wutp3nz5qxcuZLs7Gw2bdrEt99+G9oXGxtbbBxB4lmCEZEY4EXgUqA1cK2ItM5T7GZgt6qeCjwL/D1s3zpVjXdft5U0Huvk907Lei1pWa+l32EYE9KxY0f69etH+/btufLKK0lISKBOnToFlh0/fjyvv/467dq1Y8qUKTz//PMA3HLLLSxatIhOnTqxdOnSIls94AzFMmHCBC6//HK6du3KySefXGC5iy66iNWrVxMfH8+MGTPy7T///PNp0aIFZ511FqNGjcp1s8Dw4cNp164d1113XaRvha+kuKbkMVcs0hl4VFV7u+sPAKjq38LKzHPLLBGRysBvQCPgZOAjVY34Ju+EhARNTEwsdP/I176mZrVYnhjc6Zh+HlO4z9Z/BsAlLS/xORITFD/99BNnnnmmrzEcOHCAmjVrkpqayoUXXsiECRMKvLMrmhT0exGR5aqa4MX5vOyDaQxsCltPAc4trIyqZorIXqCBu6+FiKwA9gFjVPXLkgRzODOLqpWrFl/QHLW//vevgCUYEyzDhw9n9erVpKWlMXTo0KhPLn7wMsEU1NuVt7lUWJmtQDNV3SkiHYAPRKSNqu7LdbDIcGA4QLNmzYoMJs06+T0z5Y9T/A7BmHymTZvmdwhRz8tO/hSgadh6E2BLYWXcS2R1gF2qmq6qOwFUdTmwDjgt7wlUdYKqJqhqQqNGjYoM5rB18numaZ2mNK3TtPiCxpio4mWCWQa0EpEWIlIFGATMzlNmNjDUXR4AfK6qKiKN3JsEEJGWQCtgfUmCsduUvTN37Vzmrp3rdxjGmIDx7BKZ26cyApgHxAATVXWViIwFElV1NvAaMEVE1gK7cJIQwIXAWBHJBLKA21R1V0niSc/ItgctPfLkV08C0OfUPj5HYowJEk8ftFTVOcCcPNseCVtOA64u4LiZwMxSjIOMrGyq2CUyT0wfMN3vEIwxARQVf9JnZDlPxVoLxhsn1DyBE2qe4HcYxuTy/PPP07ZtW9q0acNzzz0X2r5r1y569uxJq1at6NmzJ7t37wZg5syZtGnThgsuuICdO3cCsG7dOgYNGlRg/V4pjWH4a9asWUrRlExUfONmZDoJpnJMVPy4Ze7D5A/5MPlDv8MwJiQpKYlXXnmFb7/9lu+//56PPvqINWvWAM6gkT169GDNmjX06NGDJ590LvE+88wzfPPNN9xwww2hO9DGjBnDuHHjIjpnVlaWNz9MORYV37g5LZhYa8F44pklz/DMkmf8DsOYkJ9++onzzjuP6tWrU7lyZbp168b7778POMO+DB3q3Fs0dOjQ0MjJlSpVIj09ndTUVGJjY/nyyy858cQTadWqVaHnqVmzJo888gjnnnsuS5YsYfny5XTr1o0OHTrQu3dvtm7dCsArr7xCx44dad++PVdddRWpqakA/P777/zxj3+kffv2tG/fnsWLnYn7cobhb9OmDb169QqNMbZu3Tr69OlDhw4duOCCC/j5558B2LBhA507d6Zjx448/PDDHryjx0hVK8SrQ4cOWpjtew9pr7Ef6cfL/1doGXPsth/crtsPbvc7DBMgq1evzrXe7fVu+vqK11VV9XDmYe32ejed8v0UVVU9ePigdnu9m07/cbqqqu45tEe7vd5NZ66eqarO56vb69109s+zVVV16/6tEZ2/VatWumPHDj148KCed955OmLECFVVrVOnTq6ydevWVVXVTz/9VM855xy94oordM+ePdqrVy/dtWtXkecBdMaMGc7Pdfiwdu7cWbdt26aqqtOnT9cbb7xRVVV37NgROuahhx7S8ePHq6rqNddco88++6yqqmZmZuqePXt0w4YNGhMToytWrFBV1auvvlqnTHHeq4svvlh/+eUXVVX95ptv9KKLLlJV1b59++rkyZNVVfWFF17QGjVqFPq+FPAzJKpH38tRMZpyZk4Lxi6ReaJh9YZ+h2BMLmeeeSb3338/PXv2pGbNmrRv357KlYv+uuvZsyc9e/YEYPLkyVx22WUkJyfz9NNPU69ePZ5//nmqV6+e65iYmBiuuuoqAJKTk0lKSgrVkZWVxYknngg4l+zGjBnDnj17OHDgAL179wbg888/54033gjVVadOHXbv3l3gMPwHDhxg8eLFXH31kfui0tPTAfj666+ZOdO5L2rIkCHcf//9x/7mlaKoSDA5l8gqxxzbUNqmaO/99B4AV555pc+RmKAKn8ohNiY213r12Oq51uvE1cm13rB6w1zrkd5QcvPNN3PzzTcD8OCDD9KkSRMAjj/+eLZu3cqJJ57I1q1bOe6443Idl5qayuTJk5k3bx69evVi1qxZTJs2jalTp3LLLbfkKhsXF0dMjHN3qqrSpk0blixZki+WYcOG8cEHH9C+fXsmTZrEwmKmHyhoGP7s7Gzq1q0bmnAsr2OdKsBLUfEn/ZEEExU/bpkbv3Q845eO9zsMY3LZtm0bAL/++ivvvfce1157LeBM7jV58mTAaan0798/13FPPfUUd911F7GxsRw6dAgRoVKlSqF+k8KcfvrpbN++PZRgMjIyWLVqFQD79+/nxBNPJCMjg6lTp4aO6dGjBy+99BLgtHj27duXv2JX7dq1adGiBe+88w7gJLTvv/8ecEZgDp9yICii4hvXLpF5a9agWcwaNMvvMIzJ5aqrrqJ169b07duXF198kXr16gEwevRo5s+fT6tWrZg/fz6jRx+ZqmrLli0kJiaGks59993Heeedx+TJk0OzXhamSpUqvPvuu9x///20b9+e+Pj4UKf9uHHjOPfcc+nZsydnnHFG6Jjnn3+eL774grPOOosOHTqEElJhpk6dymuvvUb79u1p06YNs2bNCtXz4osv0rFjR/bu3Xv0b5ZHPBuuv6wVNVz/qk27uHfSEp4Y3IkOpxQ9ZpkxpuSCMFy/ya+sh+uPij/pM7OcJGqXyLwxI2kGM5LyT5xkjIlu1slvSuylROca8sC2A32OxBgTJFGRYDJDQ8XYWGRemHPdnOILGWOiTlQkmFALppK1YLxQPbZ68YVM1FHVQN46G6386G+Pik6JTLtN2VNv/vAmb/7wpt9hmACJi4tj586dvnypmfxUlZ07dxIXF1em542qFozdpuyNV797FYDr213vcyQmKJo0aUJKSgrbt2/3OxTjiouLCz1sWlY8TTAi0gd4HmfCsVdV9ck8+6sCbwAdgJ3AQFXd6O57ALgZZ8Kxkao671jjyLmLzAa79Mb8IfP9DsEETGxsLC1atPA7DOMzz75x3SmPXwQuBVoD14pI6zzFbgZ2q+qpwLPA391jW+PMbtkG6AP8O2cK5WNhT/J7KzYmltiYWL/DMMYEjJffuJ2Ataq6XlUPA9OB/nnK9Acmu8vvAj3E6RXsD0xX1XRV3QCsdesrVHYR13oz7TZlT01aOYlJKyf5HYYxJmC8TDCNgU1h6ynutgLLqGomsBdoEOGxuWz4fT8TP/+ZLbsO5ttnQ8V4yxKMMaYgXvbBFNRcyNvMKKxMJMciIsOB4e5q+s09zkwqKqCqjxS1t8w0BHb4HUQEjjpOudGXFmKFfT99YnGWnvIQI8DpXlXsZYJJAZqGrTcBthRSJkVEKgN1gF0RHouqTgAmAIhIolfj6ZQmi7N0WZyly+IsPeUhRnDi9KpuL68ZLQNaiUgLEamC02k/O0+Z2cBQd3kA8Lk7w9psYJCIVBWRFkAr4FsPYzXGGFPKPGvBqGqmiIwA5uHcpjxRVVeJyFicKTpnA68BU0RkLU7LZZB77CoReRtYDWQCd6pqllexGmOMKX2ePgejqnOAOXm2PRK2nAZcnfc4d9/jwONHcboJxxKjDyzO0mVxli6Ls/SUhxjBwzgrzHwwxhhjgsXu2zXGGOOJCpFgRKSPiCSLyFoRGV38EaVyzokisk1EksK21ReR+SKyxv23nrtdRGS8G98PInJO2DFD3fJrRGRo2PYOIvKje8x4OYZhaUWkqYh8ISI/icgqEbkroHHGici3IvK9G+dj7vYWIrLUPecM92YR3Js/ZrjnXCoizcPqesDdniwivcO2l9pnRERiRGSFiHwU1DhFZKP7e1mZc5dQ0H7vbj11ReRdEfnZ/Zx2DlqcInK6+z7mvPaJyN0BjPMe9/9Pkoi8Jc7/K38/m6parl84NxCsA1oCVYDvgdZlcN4LgXOApLBtTwGj3eXRwN/d5cuAT3Ce7zkPWOpurw+sd/+t5y7Xc/d9C3R2j/kEuPQYYjwROMddrgX8gjNsT9DiFKCmuxwLLHXP/zYwyN3+H+B2d/kO4D/u8iBghrvc2v39VwVauJ+LmNL+jAD3AtOAj9z1wMUJbAQa5tkWqN+7W89k4E/uchWgbhDjDIs3BvgNODlIceI8iL4BqBb2mRzm92fT0y/hsni5v5R5YesPAA+U0bmbkzvBJAMnussnAsnu8svAtXnLAdcCL4dtf9nddiLwc9j2XOVKEO8soGeQ4wSqA98B5+I8pFY57+8Z587Ezu5yZbec5P3d55Qrzc8IzjNZC4CLgY/c8wYxzo3kTzCB+r0DtXG+FCXIceaJrRfwddDi5MjoJ/Xdz9pHQG+/P5sV4RLZUQ8r46HjVXUrgPvvce72wmIsantKAduPmdsEPhundRC4OMW57LQS2AbMx/lraY86Qwjlrftohxgqzc/Ic8D/A7Ld9QYBjVOBT0VkuTgjXkDwfu8tge3A6+JccnxVRGoEMM5wg4C33OXAxKmqm4GngV+BrTifteX4/NmsCAkmomFlfHa0Q+KU6s8kIjWBmcDdqrqvqKJHGU+pxamqWaoaj9NC6AScWUTdvsQpIlcA21R1efjmIur28/d+vqqegzOa+Z0icmERZf2KszLOZeaXVPVs4CDOpabC+P3/qArQD3inuKJHGU+J43T7f/rjXNY6CaiB87svrN4yibEiJJiIhpUpI7+LyIkA7r/b3O2FxVjU9iYFbD9qIhKLk1ymqup7QY0zh6ruARbiXLuuK84QQnnrDsUjkQ0xVFqfkfOBfiKyEWd08ItxWjRBixNV3eL+uw14HydpB+33ngKkqOpSd/1dnIQTtDhzXAp8p6q/u+tBivMSYIOqblfVDOA9oAt+fzZLcj0yCC+cv4LW42TunM6nNmV07ubk7oP5B7k7/Z5yly8nd6fft+72+jjXoOu5rw1AfXffMrdsTqffZccQn+BM6PZcnu1Bi7MRUNddrgZ8CVyB85dieAflHe7yneTuoHzbXW5D7g7K9Tidk6X+GQG6c6STP1Bx4vz1WitseTHOvEqB+r279XwJnO4uP+rGGLg43bqmAzcG8f8RTp/lKpw+TMG5eeL//P5sev4lXBYvnLs2fsG5bv9QGZ3zLZxrnRk42f1mnGuYC4A17r85Hx7BmXxtHfAjkBBWz004892szfPhTQCS3GNeIE9HaIQxdsVpxv4ArHRflwUwznbACjfOJOARd3tLnLtr1rr/Uaq62+Pc9bXu/pZhdT3kxpJM2J04pf0ZIXeCCVScbjzfu69VOfUE7ffu1hMPJLq/+w9wvniDGGd1nFl364RtC1ScwGPAz249U3CShK+fTXuS3xhjjCcqQh+MMcaYALIEY4wxxhOWYIwxxnjCEowxxhhPWIIxxhjjCUswJiqJyEIR8Xy+dBEZ6Y4SPDXP9gQRGe8udxeRLqV4zuYiMrigcxlTljyd0dKYikhEKuuR8Z2KcwfOswQbwjeqaiLO8x/gPFNzAOeByNKIoTkwGGfE57znMqbMWAvGBJb7l/hPIvKKO8/FpyJSzd0XaoGISEN3+BZEZJiIfCAiH4rIBhEZISL3uoMpfiMi9cNOcb2ILHbnz+jkHl9DnLl+lrnH9A+r9x0R+RD4tIBY73XrSRKRu91t/8F50G22iNyTp3x3EfnIHYT0NuAeceYauUBEGonITDeGZSJyvnvMoyIyQUQ+Bd5w358vReQ795XTCnoSuMCt756cc7l11Hffnx/c96NdWN0T3fd1vYiMDHs/PhZnrp4kERlYst+qiSolfaLZXvby6oXzl3gmEO+uvw1c7y4vxH1CGmgIbHSXh+E8nVwLZwiavcBt7r5ncQb8zDn+FXf5Qtwhf4Anws5RF+fJ5RpuvSm4T2vnibMDzhPbNYCaOE/Pn+3u20ieYfPd7d05MhLAo8CosH3TgK7ucjPgp7Byyzky50d1IM5dbgUk5q27gHP9C/iLu3wxsDKs7sU4T383xHlqPRa4Kud9csvVyfuz2Mtehb3sEpkJug2qutJdXo6TdIrzharuB/aLyF7gQ3f7jzjD0uR4C0BV/ysitUWkLs58H/1EZJRbJg7nSx5gvqruKuB8XYH3VfUggIi8B1yAM/zNsbgEaC1HJjWsLSK13OXZqnrIXY4FXhCReCALOC2CurviJA1U9XMRaSAiddx9H6tqOpAuItuA43Hes6dF5O84SerLY/yZTBSyBGOCLj1sOQtnMExwWjY5l3jjijgmO2w9m9yf+bzjJOUMS36VqiaH7xCRc3GGky/IMU0XXIRKOJNBHQrf6Cac8BjuAX4H2rvHpEVQd1HDrud9ryur6i8i0gFnHKq/icinqjo2op/CRD3rgzHl1UacS1MAA46xjoEAItIV2Kuqe3Fm8Ps/cb/NReTsCOr5L/AHEakuzoRZf8QZJThS+3Eu6eX4FBiRs+K2UApSB9iqqtnAEJxRbwuqL2+s17n1dgd2aBFzBInISUCqqr6JM6HVOYWVNSYvSzCmvHoauF1EFuP0GRyL3e7x/8EZDRtgHM6lpx9EJMldL5KqfgdMwhmVdinwqqoezeWxD4E/5nTyAyOBBLcjfjXOTQAF+TcwVES+wbk8ltO6+QHIdDvm78lzzKM5dePcDDC0mNjOAr4VZ7bRh4C/HsXPZaKcjaZsjDHGE9aCMcYY4wlLMMYYYzxhCcYYY4wnLMEYY4zxhCUYY4wxnrAEY4wxxhOWYIwxxnjCEowxxhhP/H9NhkPhYuMufQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "for i in range(rep):\n",
    "    plt.plot(mi_list[i,:],color='steelblue')\n",
    "    for t in range(mi_list[i].shape[0]):\n",
    "        if (mi_list[0,t]>.9*mi):\n",
    "            plt.axvline(t,label='90% reached',linestyle=':',color='green')\n",
    "            break\n",
    "#plt.title(\"Plot of MI estimates against number of iteractions\")\n",
    "plt.xlim((0,mi_list[0].shape[0]))\n",
    "plt.ylim((0,mi*1.1))\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig(fig_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
